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Abstract 

This  thesis  presents  high-order,  discontinuous  Galerkin  (DG)  discretizations  of  the  Reynolds- 
Averaged  Navier-Stokes  (RANS)  equations  and  an  output-based  error  estimation  and  mesh 
adaptation  algorithm  for  these  discretizations.  In  particular,  DG  discretizations  of  the 
RANS  equations  with  the  Spalart-Allmaras  (SA)  turbulence  model  are  examined.  The 
dual  consistency  of  multiple  DG  discretizations  of  the  RANS-SA  system  is  analyzed.  The 
approach  of  simply  weighting  gradient  dependent  source  terms  by  a  test  function  and  inte¬ 
grating  is  shown  to  be  dual  inconsistent.  A  dual  consistency  correction  for  this  discretization 
is  derived.  The  analysis  also  demonstrates  that  discretizations  based  on  the  popular  mixed 
formulation,  where  dependence  on  the  state  gradient  is  handled  by  introducing  additional 
state  variables,  are  generally  asymptotically  dual  consistent.  Numerical  results  are  pre¬ 
sented  to  confirm  the  results  of  the  analysis. 

The  output  error  estimation  and  output-based  adaptation  algorithms  used  here  are 
extensions  of  methods  previously  developed  in  the  finite  volume  and  finite  element  com¬ 
munities.  In  particular,  the  methods  are  extended  for  application  on  the  curved,  highly 
anisotropic  meshes  required  for  boundary  conforming,  high-order  RANS  simulations.  Two 
methods  for  generating  such  curved  meshes  are  demonstrated.  One  relies  on  a  user-defined 
global  mapping  of  the  physical  domain  to  a  straight  meshing  domain.  The  other  uses  a  lin¬ 
ear  elasticity  node  movement  scheme  to  add  curvature  to  an  initially  linear  mesh.  Finally, 
to  improve  the  robustness  of  the  adaptation  process,  an  “unsteady”  algorithm,  where  the 
mesh  is  adapted  at  each  time  step,  is  presented.  The  goal  of  the  unsteady  procedure  is  to 
allow  mesh  adaptation  prior  to  converging  a  steady  state  solution,  not  to  obtain  a  time- 
accurate  solution  of  an  unsteady  problem.  Thus,  an  estimate  of  the  error  due  to  spatial 
discretization  in  the  output  of  interest  averaged  over  the  current  time  step  is  developed. 
This  error  estimate  is  then  used  to  drive  an  /i-adaptation  algorithm. 

Adaptation  results  demonstrate  that  the  high-order  discretizations  are  more  efficient 
than  the  second-order  method  in  terms  of  degrees  of  freedom  required  to  achieve  a  desired 
error  tolerance.  Furthermore,  using  the  unsteady  adaptation  process,  adaptive  RANS  sim¬ 
ulations  may  be  started  from  extremely  coarse  meshes,  significantly  decreasing  the  mesh 
generation  burden  to  the  user. 
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Chapter  1 


Introduction 


In  recent  decades,  Computational  Fluid  Dynamics  (CFD)  technology  has  achieved  signifi¬ 
cant  maturity.  In  particular,  CFD  is  widely  used  throughout  industry,  academia,  and  gov¬ 
ernment  for  analysis  and  design  of  aerospace  vehicles.  Despite  this  widespread  use,  many 
challenging  problems  remain.  For  example,  when  high  accuracy  simulations  are  necessary, 
computational  costs  using  current  industry  standard  techniques  are  very  large.  Further¬ 
more,  it  is  unclear  if  the  grid  generation  practices  and  second-order  finite  volume  methods 
typically  used  by  the  aerospace  community  are  adequate  given  the  stringent  accuracy  re¬ 
quirements  of  aerodynamic  design. 

To  assess  the  current  state  of  the  art  in  CFD  for  applied  aerodynamics,  one  can  examine 
the  results  of  the  recent  American  Institute  of  Aeronautics  and  Astronautics  (AIAA)  Drag 
Prediction  Workshops  (DPW)  [70,  65,  75].  These  workshops,  held  in  June  2001,  June  2003, 
and  June  2006,  were  convened  with  the  explicit  goals  of  assessing  CFD  as  a  practical  tool  for 
computation  of  aerodynamic  forces  for  industry  relevant  geometries  and  identifying  areas 
for  additional  research  and  development. 

Over  the  three  workshops,  increasing  effort  has  been  directed  to  determining  the  effects 
of  spatial  discretization  error.  In  the  most  recent  workshop,  four  geometries  were  studied: 
two  wing-body  configurations  and  two  wing-alone  configurations  [75,  101].  The  wing-alone 
geometries  were  specifically  included  to  minimize  the  complexity  of  the  flow  physics  and 
enable  grid  convergence  studies.  The  results  show  that  the  uncertainties  associated  with 
standard  CFD  methods  are  unacceptably  high.  In  particular,  the  standard  deviation  of 
the  submitted  values  of  the  total  drag  increases  with  the  number  of  mesh  points  for  both 
wing-alone  geometries  [75].  Also,  the  magnitude  of  the  standard  deviation,  5-7  drag  counts, 
is  substantially  larger  than  the  uncertainty  desired  by  airframe  designers  [75,  100].  Further¬ 
more,  restricting  comparison  to  only  those  submissions  that  used  the  SA  turbulence  model, 
the  spread  in  the  drag  results  extrapolated  to  continuum  is  more  than  20  counts,  or  more 
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than  ten  percent  of  the  total  drag  [101]. 

Additional  evidence  that  discretization  error  plays  a  significant  role  in  the  errors  in 
current  CFD  methods  is  presented  by  Mavriplis  [73].  In  particular,  he  demonstrates  that, 
even  for  very  large  grids,  asymptotic  results  appear  to  be  different  for  different  families  of 
self-similar  grids.  This  behavior  has  been  attributed  to  the  large  range  of  scales  present 
in  high  speed,  turbulent  aerodynamic  flows.  This  large  range  of  scales  makes  it  difficult 
to  adequately  resolve  all  regions  of  the  flow  via  global  uniform  refinement  starting  from 
an  arbitrary  mesh.  Thus,  even  “fine”  meshes  may  produce  inaccurate  results  if  important 
regions  of  the  flow  are  not  well  resolved. 

These  results  show  that,  not  only  is  discretization  error  a  significant  contributor  to  the 
error  in  current  CFD  solutions,  it  can  be  very  difficult  to  detect,  even  for  expert  prac¬ 
titioners.  Thus,  there  exists  a  need  in  the  CFD  community  for  additional  research  and 
development  aimed  at  detecting  and  reducing  errors  associated  with  spatial  discretization. 
High-order,  adaptive  techniques  have  significant  promise  to  accomplish  this  aim. 

1.1  Objective 

The  objective  of  this  work  is  to  develop  a  high-order,  adaptive  method  for  the  simulation  of 
high  Reynolds  number,  turbulent  flows  and  to  demonstrate  the  performance  of  the  method 
for  two-dimensional  aerodynamic  test  cases. 

1.2  Previous  Work 

1.2.1  High-Order  Methods 

High-order  methods  have  significant  potential  to  decrease  the  impact  of  discretization  error 
on  the  accuracy  of  CFD  solutions.  Most  CFD  methods  in  widespread  use  in  the  aerospace 
industry  achieve,  at  best,  E  cc  where  FI  is  a  measure  of  the  error  in  the  solution 
and  h  is  a  measure  of  the  mesh  spacing.  Thus,  in  the  context  of  this  work,  high-order 
methods  are  those  that  achieve  E  (x  K' ,  where  r  >  2,  for  sufficiently  smooth  problems. 
While  second-order  finite  volume  discretizations  are  popular,  high-order  methods  have  been 
extensively  studied.  These  efforts  have  led  to  many  types  of  high-order  schemes,  including 
finite  difference,  finite  volume,  and  finite  element  methods. 

In  the  context  of  finite  difference  methods,  Lele  [68]  introduced  up  to  tenth-order 
compact  finite  difference  schemes.  This  work  was  extended  and  applied  by  Visbal  and 
Gaitonde  [104],  who  used  high-order  compact  difference  methods  to  solve  the  compress¬ 
ible  Navier-Stokes  equations  on  curvilinear  meshes.  Additional  work  in  high-order  finite 
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difference  methods  for  aerodynamics  was  conducted  by  Zingg  et  al.  [113],  who  showed 
that  high-order  compact  difference  methods  are  more  efficient,  in  terms  of  number  of  nodes 
required  to  accurately  compute  drag,  than  typical  second-order  finite  differences. 

While  these  high-order  finite  difference  techniques  have  been  successfully  used  in  many 
cases,  their  application  is  limited  to  structured  meshes.  To  minimize  mesh  generation  ef¬ 
fort  for  complex  geometries,  unstructured  meshes  are  of  interest.  For  unstructured  meshes, 
Barth  [11]  pioneered  high-order  finite  volume  methods  using  the  least-squares  /c-exact  re¬ 
construction  method.  For  this  method,  additional  nodes  and  control  volumes  are  added 
to  each  element  to  enable  high-order  reconstruction  without  extending  the  reconstruction 
support  outside  of  the  element.  More  recently,  a  similar  idea,  known  as  the  spectral  volume 
method,  has  been  developed  by  Wang  [105].  In  this  scheme,  each  mesh  cell — i.e.  spectral 
volume — is  divided  into  sub-cells.  The  state  averages  on  these  sub-cells  are  then  used  to 
build  a  high-order  reconstruction  of  the  solution  within  the  spectral  volume. 

In  the  context  of  finite  element  methods,  researchers  in  the  early  1980s  pioneered  the 
so-called  p-type  finite  element  method.  In  the  p-type  method,  the  grid  spacing,  h,  remains 
fixed  while  the  interpolation  order,  p,  is  increased  to  improve  the  accuracy  of  the  solution. 
Babuska  et  al.  [8]  applied  the  p-type  method  to  the  elasticity  equations  in  1981  and  con¬ 
cluded  that  the  p-type  method  achieved  superior  performance  in  terms  of  the  degrees  of 
freedom  required  to  achieve  a  desired  accuracy.  Later,  Patera  [81,  64]  introduced  a  vari¬ 
ant  of  the  p-type  method,  known  as  the  spectral  element  method,  and  used  it  to  solve  the 
incompressible  Navier-Stokes  equations. 

While  finite  element  methods  offer  a  conceptually  simple  path  to  high-order  accuracy, 
it  is  well  known  that  the  standard,  continuous  Galerkin  method  is  inappropriate  for  use 
on  convection-dominated  problems.  This  drawback  stems  from  the  fact  that  the  basic  con¬ 
tinuous  Galerkin  discretization  is  unstable  for  convection.  Thus,  even  for  subsonic  flows, 
the  discretization  of  the  Euler  or  Navier-Stokes  equations  with  finite  element  methods  re¬ 
quires  the  addition  of  stabilization.  One  popular  technique  is  the  Streamline-Upwind  Petrov 
Galerkin  method  [61].  Another  method  that  has  received  significant  attention  is  the  dis¬ 
continuous  Galerkin  (DG)  method. 

1.2.2  Discontinuous  Galerkin  Methods 

For  achieving  high-order  accuracy,  DG  is  attractive  for  two  reasons.  First,  it  allows  the  de¬ 
velopment  of  stable,  high-order  accurate  discretizations  of  convection-dominated  problems 
using  upwinding  methods  developed  in  the  finite  difference  and  finite  volume  communities. 
Second,  the  resulting  discretizations  have  element-wise  compact  stencils.  These  compact 
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stencils  simplify  the  task  of  achieving  high-order  accuracy  on  unstructured  meshes  and  near 
boundaries,  and  they  allow  the  development  of  efficient  solution  methods. 

In  1973,  Reed  and  Hill  [90]  introduced  the  DG  method  for  the  neutron  transport  equa¬ 
tion.  Since  that  time,  there  has  been  a  rapid  proliferation  of  new  DG  techniques  as  well  as 
analyses  and  applications  of  those  techniques.  Gockburn  et  al.  [26]  provides  an  excellent 
review  of  work  in  DG  methods  through  the  year  2000.  Highlights  of  this  review  as  well  as 
relevant  recent  advances  are  summarized  here. 

In  1974,  LeSaint  and  Raviart  [69]  proved  the  first  a  priori  error  estimates  of  the  DG 
method  for  linear  hyperbolic  problems.  They  showed  that,  assuming  a  smooth  exact  so¬ 
lution,  the  error  measured  in  the  L^-norm  is  0{h^).  Later,  Johnson  and  Pitkaranta  [63] 
and  Richter  [92]  improved  upon  the  estimate  of  LeSaint  and  Raviart.  Specifically,  Johnson 
and  Pitkaranta  proved  that,  in  the  most  general  case,  is  the  optimal  convergence 

rate  in  L^,  while  Richter  showed  that,  assuming  the  characteristic  direction  is  not  exactly 
aligned  with  the  grid,  0{h^^^)  can  be  obtained. 

The  extension  from  the  method  of  Reed  and  Hill  for  linear  problems  to  nonlinear  hyper¬ 
bolic  problems  is  accomplished  by  the  use  of  a  Riemann  solver  to  evaluate  the  flux  across 
element  boundaries.  Riemann  solvers  have  been  extensively  developed  in  the  finite  volume 
community  [93,  97].  The  first  application  of  the  DG  method  to  a  nonlinear  hyperbolic 
problem  was  accomplished  by  Ghavent  and  Salzano  [24]  using  Godunov’s  flux. 

A  breakthrough  in  the  application  of  DG  methods  to  nonlinear  hyperbolic  problems 
using  explicit  time  integration  was  made  by  Gockburn,  Shu,  and  co-authors  [28,  27,  25,  30], 
who  introduced  the  Runge  Kutta  Discontinuous  Galerkin  (RKDG)  method.  The  original 
RKDG  method  uses  an  explicit  TVD  second-order  Runge  Kutta  scheme  introduced  by  Shu 
and  Osher  [94].  This  method  was  later  generalized  to  be  high-order  accurate  in  time  as  well 
as  space. 

Independent  of  the  above  work,  Allmaras  [2]  and  Allmaras  and  Giles  [4]  developed  a 
second-order  DG  scheme  for  the  2-D  Euler  equations.  Their  method  is  the  extension  of 
van  Leer’s  method  of  moments  [98]  from  the  1-D,  linear  wave  equation  to  the  2-D  Euler 
equations.  Thus,  it  requires  that  state  and  gradient  averages  be  computed  at  each  cell  to 
allow  linear  reconstruction  of  the  state  variables.  Halt  [49]  later  extended  this  technique  to 
be  higher-order  accurate. 

For  elliptic  operators,  in  the  late  1970s  and  early  1980s,  Arnold  [6]  and  Wheeler  [107] 
introduced  discontinuous  finite  element  methods  known  as  penalty  methods.  More  recently, 
many  researchers  [14,  78,  29,  15,  16,  31,  20]  have  applied  DG  methods  to  diffusion  problems. 
One  procedure,  pioneered  by  Bassi  and  Rebay  [14,  15]  and  generalized  by  Gockburn  and 
Shu  [29,  31],  is  to  rewrite  a  second-order  equation  as  a  first-order  system  and  then  discretize 
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the  first-order  system  using  the  DG  formulation.  Methods  derived  in  this  fashion  will 
be  referred  to  as  mixed  formulations.  The  first  mixed  formulation  developed  by  Bassi 
and  Rebay  (BRl)  is  not  coercive  and  produces  an  extended  stencil.  However,  a  slight 
modification  of  BRl  leads  to  a  coercive  scheme  with  a  nearest-neighbor  stencil,  referred  to  as 
the  second  method  of  Bassi  and  Rebay  (BR2).  The  generalization  of  the  first-order  system 
idea  by  Cockburn  and  Shu  leads  to  the  so-called  Local  Discontinuous  Galerkin  methods 
(LDG).  In  multiple  dimensions,  LDG  has  an  extended  stencil,  but  a  recent  modification 
of  LDG  by  Persson  and  Peraire  [82]  known  as  Gompact  Discontinuous  Galerkin  (CDG) 
recovers  a  compact  stencil  while  retaining  the  attractive  properties  of  LDG. 

The  penalty-type  and  mixed  formulation  DG  methods  have  been  brought  into  a  single 
analysis  framework.  This  framework,  introduced  by  Arnold  et  al.  [7],  provides  for  a  unified 
analysis,  including  error  estimates,  of  these  schemes.  One  recently  introduced  scheme  which 
has  not  been  incorporated  into  this  framework  is  that  of  van  Leer  [99],  who  uses  a  patch 
reconstruction  of  the  solution  to  evaluate  the  flux  across  element  boundaries. 

This  work  examines  three  DG  discretizations  of  the  RANS  equations  with  the  Spalart- 
Allmaras  (SA)  turbulence  model.  The  application  of  DG  discretizations  to  the  RANS 
equations  has  been  somewhat  limited.  At  the  time  of  this  writing,  the  author  is  aware 
of  three  DG  implementations  of  the  RANS  equations.  Specifically,  Bassi  and  Rebay  [12] 
used  the  BR2  method  to  discretize  the  RANS  equations  coupled  with  the  k  —  uo  turbulence 
model,  and  Nguyen,  Persson,  and  Peraire  [76]  used  GDG  for  the  RANS  equations  coupled 
with  the  SA  turbulence  model.  Most  recently,  Landmann  [66]  has  applied  both  LDG  and 
BR2  to  discretize  the  RANS  equations  coupled  with  the  SA  and  k  —  oj  turbulence  models. 
All  three  of  these  methods  are  based  on  mixed  formulations,  which  are  examined  in  detail 
in  Chapter  3. 

1.2.3  Error  Estimation  and  Adaptation 

It  is  widely  recognized  that  error  estimation  and  adaptation  increase  the  usefulness  of  CFD 
computations.  Estimating  the  error  in  a  CFD  solution  and  generating  appropriate  high- 
quality  meshes  are  difficult  tasks,  even  for  experts  in  CFD  and  aerodynamics.  By  enabling 
the  user  to  set  the  error  tolerance  of  the  computation  and  automatically  adapting  the 
discretization  to  satisfy  that  tolerance,  the  confidence  in  the  computed  solution  increases 
and  the  user  is  freed  from  the  onerous  task  of  generating  a  high-quality  mesh. 

Given  the  potential  impact  of  error  estimation  and  adaptation  on  the  usefulness  of 
CFD,  many  researchers  have  studied  error  estimation  and  adaptation  algorithms.  A  brief 
summary  of  relevant  advances  is  given  here. 
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Error  Estimation 


This  section  gives  a  brief  review  of  error  estimation  techniques  used  to  drive  mesh  adap¬ 
tation  algorithms.  Thus,  the  techniques  discussed  include  both  rigorous  a  posteriori  error 
estimation  algorithms  as  well  as  more  ad  hoc  ideas  used  to  construct  “error  indicators”  to 
target  mesh  adaptation. 

A  common  technique  for  driving  adaptation  in  the  CFD  community  is  feature  detection. 
This  method  aims  to  identify  the  dominant  flow  features,  such  as  shock  waves  and  boundary 
layers,  by  finding  regions  with  large  gradients  [9,  86,  106].  The  underlying  assumption  is  that 
the  error  in  these  regions  dominates  the  error  in  the  solution.  Thus,  the  identified  large 
gradient  regions  are  targeted  for  grid  refinement.  Such  methods  have  been  successfully 
applied  to  some  flow  problems,  but  they  are  clearly  ad  hoc.  Furthermore,  it  has  been 
demonstrated  that  simply  refining  the  dominant  features  of  the  flow  can  lead  to  incorrect 
results  [106].  For  example,  small  errors  in  the  upstream  flow  can  affect  shock  or  separation 
locations,  leading  to  large  errors  in  computed  outputs.  These  errors  are  not  reduced  by 
continually  refining  the  shock  or  boundary  layer. 

A  less  ad  hoc  procedure  is  offered  by  Zienkiewicz  and  Zhu  [110,  111,  112],  who  propose 
an  error  estimation  algorithm  based  on  recovery.  The  idea  underlying  these  techniques 
is  to  use  the  current  discrete  solution  to  reconstruct  a  better  approximation  of  the  exact 
solution.  Then,  the  difference  between  the  discrete  solution  and  the  reconstruction  can  be 
used  to  assess  the  local  error  in  the  solution.  Many  researchers  have  used  similar  schemes 
based  on  interpolation  error  estimates  [21,  48,  83,  106].  However,  such  schemes  have  many 
drawbacks.  For  example,  Ainsworth  and  Oden  [1]  present  a  second-order  ODE  case  where 
the  recovery-based  error  estimate  is  zero  while  the  actual  error  can  be  arbitrarily  large. 
Furthermore,  given  the  local  nature  of  such  error  estimates,  they  may  fail  to  correctly 
capture  propagation  of  error  for  convection-dominated  problems.  This  failure  can  lead  to 
similar  problems  as  those  described  for  purely  feature-based  algorithms. 

Another  type  of  error  estimation  technique  relies  on  the  residual.  These  estimates  are 
motivated  by  the  observation  that  it  is  often  possible  to  bound  the  error  by  an  appropriately 
defined  residual  norm.  For  example,  consider  a  linear  PDE  of  the  form 

Cu  =  f. 

Then,  given  an  approximate  solution,  Uh,  the  error,  e/^  =  u  —  tt/i,  is  governed  by 


^eh  =  rh 
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where  rh  ^  f  —  Cuh-  Furthermore,  this  provides  an  error  bound  of  the  form 

||e/i||i  <  C\\rh\\2, 

where  ||  •  ||i  and  ||  •  II2  are  appropriate  norms,  and  (7  is  a  constant. 

Certainly,  the  analysis  is  more  difficult  for  nonlinear  problems,  but  in  many  cases,  it 
is  possible  to  bound  the  error  using  an  appropriate  residual  and  norm  [109,  87].  Thus, 
it  is  possible  to  construct  error  indicators  based  on  the  residual.  An  interesting  study 
is  provided  by  Zhang  et  al.  [109].  They  find  that,  for  driving  mesh  adaptation  in  one¬ 
dimensional  subsonic  flow,  both  recovery-based  and  residual-based  error  estimation  schemes 
are  adequate,  but  that  the  residual-based  method  is  more  efficient  for  fine  meshes.  For  one¬ 
dimensional  transonic  flow  with  shocks,  they  conclude  that  the  residual-based  method  is 
superior  because  it  is  able  to  account  for  the  transport  of  error  due  to  convection.  However, 
neither  estimate  is  adequate  for  driving  adaptation  in  two  dimensions. 

A  significant  drawback  of  both  recovery-  and  residual-based  procedures  discussed  thus 
far  is  that  they  aim  to  estimate  some  global  solution  error — e.g.  the  or  norm  of 
the  error  over  the  entire  domain.  Alternatively,  the  most  important  errors  are  generally 
those  in  quantities  of  engineering  interest — e.g.  lift,  drag,  etc.  Hence,  another  class  of 
methods  has  been  developed  to  estimate  the  error  in  the  output  of  interest  directly.  Such 
algorithms  are  known  as  output-based  or  goal-oriented  error  estimation  schemes.  These 
schemes  are  generally  divided  into  two  types;  Type  I  methods,  which  depend  explicitly  on 
the  solution  of  an  appropriate  dual  problem,  and  Type  H  methods,  where  the  dependence  on 
the  dual  problem  is  eliminated.  While  both  types  of  methods  can  be  used  to  construct  error 
indicators  to  target  adaptation,  Type  I  methods  have  been  shown  to  be  superior  [55,  54] 
for  efficient  and  accurate  computation  of  functional  outputs.  Thus,  Type  I  methods  are 
examined  here.  In  particular,  the  method  used  in  this  work  is  based  on  the  Dual  Weighted 
Residual  (DWR)  method  due  to  Becker  and  Rannacher  [17,  18].  The  DWR  method  uses  the 
property  of  Galerkin  orthogonality  of  finite  element  discretizations,  combined  with  duality 
concepts,  to  express  the  output  error  in  terms  of  weighted  residuals. 

Many  implementations,  variations,  and  extensions  of  the  DWR  methods  appear  in  the 
literature.  For  example.  Pierce  and  Giles  [85,  43,  45]  and  Venditti  and  Darmofal  [102,  103] 
have  developed  duality  based  error  correction  methods  that  extend  the  super  convergence 
properties  of  finite  element  methods  to  more  general  discretizations.  Venditti  and  Dar¬ 
mofal  [103]  also  used  an  estimate  of  the  error  remaining  after  correction  to  drive  a  mesh 
adaptation  procedure.  This  procedure  was  the  first  output-based,  anisotropic  adaptation 
method  applied  to  the  RANS  equations.  The  results  show  that,  for  a  second-order  finite 
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volume  discretization,  the  output-based  approach  is  significantly  more  reliable  than  a  purely 
feature-based  adaptation  scheme. 

For  DG  finite  element  methods,  Hartmann  and  Houston  [55,  56,  57],  Hartmann  [52], 
Houston  and  Siili  [59,  60],  Lu  [72],  and  Fidkowski  [42,  41]  have  all  used  techniques  based 
directly  on  or  fundamentally  similar  to  the  DWR  method.  Similar  ideas  are  also  used  in 
this  work.  They  are  presented  in  detail  in  Chapter  4. 

Adaptation 

Given  an  error  indicator  or  estimate,  the  goal  of  the  adaptation  procedure  is  to  modify 
the  discretization  to  decrease  the  error.  For  finite  element  methods,  this  can  be  done  in 
one  of  three  ways:  p-adaptation,  where  the  order  of  the  elements  is  modified  on  a  constant 
mesh;  h-adaptation,  where  the  mesh  is  changed  but  the  order  of  the  elements  remains 
hxed;  or  /ip-adaptation,  where  both  h  and  p  are  allowed  to  change.  While  p-adaptation 
is  known  to  be  more  efficient  for  sufficiently  regular  solutions,  /i-adaptation  easily  allows 
for  the  generation  of  highly  anisotropic  meshes,  which  are  critical  for  efficiency  in  high  Re 
flows.  While  hp-adaptation  might  be  the  most  efficient,  in  that  case,  one  must  contend  with 
the  additional  difficulty  of  deciding  between  h  and  p  adaptation.  Some  researchers  have 
proposed  algorithms  for  making  this  choice  [60],  but  it  is  not  a  solved  problem.  Thus,  for 
simplicity,  /i-adaptation  alone  is  chosen  for  this  work. 

The  h-adaptation  is  driven  by  an  output-based  error  estimate.  However,  as  noted  earlier, 
appropriately  anisotropic  meshes  are  crucial  for  efficient  simulation  of  high  Re  flows.  In  this 
work,  the  desired  anisotropy  calculation  is  based  on  equidistributing  interpolation  error. 
This  approach  is  similar  to  that  of  Peraire  et  al.  [83],  who  use  the  Hessian  of  the  density 
to  determine  the  mesh  spacing  request.  Gastro-Diaz  et  al.  [21]  and  Habashi  [48]  also 
use  the  Hessian  of  a  scalar  quantity  to  determine  the  desired  mesh  spacing.  Alternatively, 
Venditti  [102,  103]  uses  an  output  error  estimate  to  determine  the  absolute  mesh  size  request 
and  only  uses  the  Hessian  to  determine  the  desired  anisotropy.  Fidkowski  [42,  41]  extended 
Venditti’s  approach  to  higher-order  discretization  methods.  A  similar  technique  based  on 
higher-order  derivatives  has  been  applied  by  Leicht  and  Hartmann  [67],  who  also  consider 
an  anisotropic  adaptation  indicator  based  on  the  inter-element  jumps  in  the  DG  solution. 
The  method  of  Fidkowski,  described  in  Ghapter  5,  is  used  here. 

1.3  Thesis  Overview 

This  thesis  describes  the  development  of  a  high-order,  h-adaptive,  DG  discretization  for  the 
RANS  equations  coupled  with  the  SA  turbulence  model.  In  particular,  the  thesis  makes 
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the  following  contributions; 


•  an  analysis  of  the  dual  consistency  of  DG  discretizations  of  source  terms  depending 
on  state  gradients,  like  those  appearing  in  the  SA  model, 

•  extension  of  the  output-based  error  estimation  implementation  to  curved  meshes, 

•  development  of  curved  mesh  generation  techniques,  including  a  global  mapping  ap¬ 
proach  and  a  linear  elasticity  mesh  movement  approach, 

•  application  of  the  output-based  error  estimation  and  adaptation  algorithms  to  high- 
order  discretizations  of  the  RANS  equations  on  curved  meshes,  and 

•  development  of  an  unsteady  adaptation  algorithm  to  improve  the  robustness  of  adap¬ 
tation  for  steady  state  problems. 

The  thesis  begins  with  a  review  of  the  RANS  equations  and  the  SA  turbulence  model 
in  Chapter  2.  Chapter  2  also  describes  modifications  to  the  SA  model  made  to  improve 
robustness.  Test  results  demonstrate  that,  as  the  mesh  is  refined,  the  modifications  to  the 
model  have  a  diminishing  effect  on  the  computed  solution.  Finally,  the  chapter  concludes 
by  detailing  the  discretizations  of  the  RANS-SA  system  used  in  this  work. 

An  analysis  of  the  dual  consistency  of  these  discretizations  is  given  in  Chapter  3.  The 
analysis  shows  that,  for  source  terms  depending  on  the  gradient  of  the  state,  the  straight¬ 
forward  method  of  weighting  the  source  term  by  a  test  function  and  integrating  leads  to  a 
dual  inconsistent  scheme.  A  dual  consistency  correction  to  the  standard  weighting  scheme 
is  derived.  Further  analysis  demonstrates  that  discretizations  based  on  the  popular  mixed 
formulation  are,  in  general,  asymptotically  dual  consistent.  Model  problem  and  RANS 
results  confirm  the  conclusions  of  the  analysis. 

Chapter  4  details  the  error  estimation  algorithm,  including  implementation  modifica¬ 
tions  used  here  to  improve  the  performance  of  the  algorithm  on  curved  meshes.  Numerical 
results  show  that  the  error  estimate  is  quite  accurate  for  two  high  Re,  two-dimensional  test 
flows.  Chapter  5  discusses  the  adaptation  algorithm,  including  two  curved  mesh  generation 
techniques.  The  first  of  these  techniques  relies  on  a  global  mapping  from  a  straight  mesh¬ 
ing  space  to  the  curved  physical  domain.  While  this  method  is  successfully  demonstrated 
in  two  dimensions,  it  is  quite  restrictive.  Thus,  a  more  general  technique,  using  a  mesh 
movement  algorithm  based  on  a  linear  elasticity  analogy,  is  developed.  This  technique  is 
demonstrated  on  two-dimensional  test  problems,  but  it  can  be  extended  to  three-dimensions 
in  a  straightforward  manner.  Chapter  5  concludes  with  multiple  two-dimensional  test  cases 
demonstrating  the  performance  of  the  adaptive  algorithm. 
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Chapter  6  provides  an  unsteady  adaptation  algorithm  for  improving  the  robustness  of  the 
algorithm  from  Chapter  5.  A  major  weakness  of  the  standard  algorithm  is  the  need  to  obtain 
a  steady  state  solution  prior  to  adaptation.  Especially  when  starting  from  coarse  meshes, 
this  need  places  an  unacceptably  high  robustness  requirement  on  the  flow  solver.  Thus, 
the  goal  of  the  new  adaptation  approach  is  to  allow  mesh  adaptation  prior  to  obtaining  a 
converged  steady  state  solution.  The  new  approach  is  referred  to  as  the  unsteady  algorithm 
because  the  mesh  is  adapted  at  every  time  step  while  marching  to  the  steady  state  solution. 
However,  the  goal  is  not  a  time  accurate  solution.  Simple  test  cases  demonstrate  that  the 
unsteady  algorithm  is  capable  of  obtaining  accurate  steady  state  RANS  solutions  starting 
from  very  coarse,  inviscid-style  meshes. 

The  thesis  ends  with  Chapter  7,  which  gives  conclusions  and  suggestions  for  future  work. 


28 


Chapter  2 


Discretization  of  the  RANS 
Equations 


This  chapter  describes  the  RANS  equations,  the  SA  turbulence  model,  and  the  discretization 
of  the  RANS-SA  system.  R  begins  with  a  brief  review  of  the  RANS  equations  and  the  SA 
turbulence  model  in  Sections  2.1  and  2.2.  Sections  2.3  and  2.4  show  the  spatial  and  temporal 
discretizations  used  in  this  work.  An  in-depth  analysis  of  the  spatial  discretizations  is  given 
in  Chapter  3.  Finally,  Section  2.5  provides  a  brief  overview  of  the  solution  method  used  to 
solve  the  discrete  system. 


2.1  The  RANS  Equations 


The  RANS  equations  are  derived  by  averaging  the  Navier-Stokes  equations.  Specifically, 
for  compressible  flows,  the  Favre  averaging  procedure  is  used.  This  procedure  as  well  as 
simplifying  assumptions  used  here  are  described  in  Appendix  A.  The  form  of  the  RANS 
equations  used  in  this  work  is  given  by 
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where  p  denotes  the  density,  Ui  are  the  velocity  components,  p  is  the  pressure,  e  is  internal 
energy,  h  is  the  enthalpy,  T  is  the  temperature,  Sij  =  ^  is  the  strain-rate  tensor, 

p  is  the  dynamic  viscosity,  pt  is  the  dynamic  eddy  viscosity,  Pr  is  the  Prandtl  number,  and 
Prt  is  the  turbulent  Prandtl  number.  The  (•)  and  (•)  notation  indicates  Reynolds-averaging 
and  Favre-averaging,  respectively.  These  averages  are  defined  in  Appendix  A. 

Equations  (2.1)  through  (2.3)  contain  more  unknowns  than  equations.  In  particular, 
the  eddy  viscosity,  pt,  which  relates  the  mean  flow  viscous  stresses  to  the  apparent  stresses 
due  to  turbulent  fluctuations,  cannot  yet  be  computed.  Thus,  to  close  the  system,  (2.1) 
through  (2.3)  are  augmented  by  the  SA  turbulence  model,  described  in  Section  2.2. 

To  avoid  confusion,  the  (•)  and  (•)  notation  is  dropped  for  the  remainder  of  the  thesis. 
All  uses  of  the  standard  flow  variables  refer  to  the  appropriate  mean  flow  quantities — e.g. 
p  is  the  Reynolds-average  density  and  Ui  is  the  Favre-average  velocity. 


2.2  The  SA  Turbulence  Model 

The  closure  of  the  RANS  system  is  accomplished  by  the  addition  of  a  turbulence  model.  In 
this  work,  the  SA  turbulence  model  [95]  is  used.  The  model  is  widely  used  in  the  aerospace 
industry  and  is  generally  regarded  as  robust.  Moreover,  it  has  been  shown  to  be  accurate 
for  most  attached  and  mildly  separated  aerodynamic  flows  [46,  89,  108,  23]. 


2.2.1  Baseline  Model 

The  model  takes  the  form  of  a  PDE  for  a  working  variable,  P,  which  is  algebraically  related 
to  the  eddy  viscosity,  pt-  In  particular,  the  eddy  viscosity  is  given  by 

l^t  =  pi>fvi, 
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where 
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S  is  the  magnitude  of  the  vorticity,  and 
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The  remaining  closure  functions  are 
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where  d  is  the  distance  to  the  nearest  wall,  Cbi  =  0.1355,  a  =  2/3,  Cb2  =  0.622,  k  =  0.41, 
Cwi  =  Cbi/K^  +  il  +  Cb2)/cr,  Cw2  =  0.3,  Cw3  =  2,  Cyi  =  7.1,  Cy2  =  0.7,  Cy3  =  0.9,  and  Pn  =  0.9. 


Some  clarifying  remarks  are  in  order.  To  begin,  the  laminar  suppression  and  trip  terms 
from  the  original  model  are  omitted  because  they  are  not  used  in  this  work.  All  cases  are 
run  fully  turbulent  with  no  effort  to  model  transition  or  force  the  flow  to  transition  at  a 
desired  location. 


In  addition,  the  form  of  the  SA  model  shown  in  (2.4)  has  been  modified  in  two  ways 
from  that  given  in  [95].  First,  it  is  a  straightforward  generalization  of  the  original  model 
to  the  compressible  flow  case.  For  incompressible  flows,  the  two  models  are  exactly  the 
same.  No  effort  was  made  in  this  work  to  optimize  the  form  of  the  model  for  flows  where 
compressibility  effects  are  highly  important.  Thus,  it  is  likely  that  another  form  would  be 
more  appropriate  for  such  cases.  For  example,  Catris  and  Aupoix  [22]  propose  a  form  of  the 
model  that  is  compatible  with  density  variations  in  the  log  layer  of  a  compressible  boundary 
layer. 

Second,  the  form  of  the  production  term  has  been  changed.  In  the  original  production 
term,  S  is  given  by  simply  S  =  S  +  S.  This  form  can  cause  robustness  problems  because 
the  production  can  be  negative,  even  for  positive  values  of  u.  Negative  values  can  occur 
because  the  fy2  closure  function  is  negative  for  approximately  1.00  <  x  <  18.40,  as  shown 
in  Figure  2-1. 

The  new  form  of  the  production  was  developed  by  Johnson  and  Allmaras  to  avoid  the 
possibility  of  negative  S  for  positive  i)  [3].  In  fact,  S  is  non-negative  for  all  f',  and  as 
S  — oo,  S  ^  {1  —  Cy3)S.  Furthermore,  the  function  is  continuous,  with  the  value  and 
derivative  matching  at  5  =  —0^28,  where  5  =  (1  —  0^2)8. 
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Figure  2-1:  SA  model  fv2  closure  function  versus  x 


2.2.2  Negative  u  Modifications 


The  exact  solution  of  (2.4)  is  non-negative,  which  agrees  with  the  physical  intuition  that  nt  > 
0.  However,  the  discrete  solution  of  (2.4)  may  not  share  this  property.  More  importantly, 
negative  D  values  can  adversely  impact  the  iterative  convergence  of  the  discrete  solution, 
even  causing  it  to  diverge  in  some  cases. 

To  ameliorate  this  behavior,  changes  can  be  made  to  the  model  for  negative  9  values. 
The  most  obvious  change  is  to  the  definition  of  the  eddy  viscosity.  The  modified  definition 
is  given  by 


pi^fvi  9  >  0 
0  9<0. 


(2.5) 


Clearly,  this  modification  ensures  that  the  eddy  viscosity  is  non-negative.  Furthermore,  the 
definition  is  continuous  and  has  continuous  first  and  second  derivatives. 

The  remaining  changes  are  motivated  by  examining  the  energy  of  the  turbulence  model 
working  variable,  One  can  derive  a  governing  equation  for  by  multiplying  (2.4) 

by  9.  The  resulting  equation  is 
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where  rj  =  fj,  +  p9,  P  =  ChiSp9,  D  =  c^ifwPP'^ /d‘^ ■ 

For  9  <  0,  the  right-hand  side  of  (2.6)  should  act  to  dissipate  et>.  This  property  ensures 
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that,  in  regions  where  9  <  0,  the  energy  of  the  turbulence  model  working  variable  will 
decrease  with  time. 


To  make  this  notion  more  precise,  define  the  integrated  energy  by 


Jn 


where  the  bounded  set  C  is  the  domain  of  interest.  Furthermore,  define  the  following 
sub-domains; 


=  {x  G  n  I  i>(x,  t)  >  0}, 
=  {x  G  n  I  i>(x,  f)  <  0}. 


Then,  E^,  =  E'd  -|-  E-  ,  where 
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To  bound  the  solution  in  regions  where  P  <  0,  consider  the  integrated  energy  contained 
in  E^ .  Given  an  initial  value  at  time  to,  one  can  determine  if  this  energy  will  grow  by 
examining  the  derivative  In  particular. 


dE^ 
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where  v  denotes  the  velocity  of  the  movement  of  the  boundary  of  and  n  is  the  outward 
pointing  unit  normal  vector. 


Consider  the  case  where  Ddil.  =  0.  Then,  assuming  that  9  is  continuous,  it  is  clear 
that  =  0,  which  implies  that  eplg^^-  =  0.  Thus, 


pe^v  ■  fids  =  0. 
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Furthermore,  using  (2.6) 
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Applying  the  divergence  theorem  and  recalling  that  =  0  gives 
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If  the  right-hand  side  of  (2.7)  is  negative  or  zero,  then  is  a  non-increasing  function  of 
time.  Thus,  E'^  is  bounded  by  its  initial  value,  which  provides  a  bound  on  the  size  of  the 
solution  in  regions  where  P  <  0. 


However,  for  the  baseline  model,  can  be  positive.  The  contribution  of  the  term 

^  {cb2pi^  —  p)  ■§^-§^  is  positive  whenever  {cb2pP  —  p)  is  positive,  which  occurs  for  y  < 
l/(cf,2  —  1)  ~  —2.65.  Furthermore,  the  contribution  of  the  production  term,  given  by 

PP  =  CbiSpp'^, 

is  non-negative  regardless  of  P.  Finally,  the  contribution  of  the  destruction  term  is 

_~n  -  _  f  EEI 

—  Cwljw  ^2  ' 

The  sign  of  this  term  is  not  entirely  determined  by  P  because  the  function  changes  sign 
at  r  ss  —1.18,  and  r  depends  on  P,  S,  and  d.  However,  for  r  <  —1.18,  is  positive,  which, 
for  P  <  0,  implies  that  —PD  is  positive. 


To  fix  these  problems,  one  can  change  the  model  when  P  <  0  to  ensure  that 

{cb2pP  -p)  <  0,  (2.8) 

P{P-D)  <  0.  (2.9) 

To  satisfy  these  properties,  changes  to  the  diffusion  coefficient,  p,  as  well  as  the  produc¬ 
tion,  P,  and  destruction,  D,  terms  have  been  made.  These  changes  were  suggested  by 
Allmaras  [3]. 
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To  satisfy  (2.8)  the  diffusion  coefficient  is  re-defined  as 


V  = 


h(i  +  x),  X>0 

+  x  +  x<o. 


(2.10) 


The  diffusion  coefficient  and  the  quantity  Cb2X  —  &re  plotted  in  Figure  2-2.  Clearly, 
the  modification  ensures  that  the  diffusion  coefficient  is  always  positive  and  Cf,2X  “  is 
always  negative.  Also,  rj  is  continuous  and  has  continuous  first  derivatives. 

The  modified  production  is  given  by 


P 


CbiSpi>,  X  >  0 

CbiSpugn,  X  <  0, 


(2.11) 


where 

_  ,  lOOOx^ 

i/n  —  1  ..  2  • 

1  +  X^ 

The  function  gn  is  chosen  such  that  P  has  continuous  first  derivatives  at  P  =  0.  This 
choice  implies  that  gn  >  0,  and  thus  vP  >  0,  for  some  negative  v.  However,  gn  <  ^  for 
X  <  — -\/l/999  ss  —0.032,  which  implies  that  vP  <  0  for  x  <  —0.032. 

Finally,  the  modified  destruction  is  given  by 


D 


c 

^WlJ  W~^ 


-r 


X>0 
X  <  0. 


(2.12) 


Clearly,  D  is  negative  for  F  <  0.  Thus,  —vD  <  0  for  F  <  0.  Furthermore,  D  is  continuous 
and  has  continuous  first  derivatives. 


2.2.3  Comparison  of  SA  Model  Versions 

To  demonstrate  the  effects  of  the  modifications  to  the  SA  turbulence  model  on  computed 
flow  solutions,  an  example  case  is  solved  using  three  versions  of  the  model:  the  original 
model  as  presented  in  [95],  the  baseline  model  presented  in  Section  2.2.1,  and  the  modified 
model  presented  in  Section  2.2.2.  For  all  the  models,  the  non- negative  form  of  given 
in  (2.5),  is  used,  as  this  modification  was  required  to  obtain  converged  solutions  for  the 
original  and  baseline  models. 

The  case  is  Moo  =  0.25,  Rec  =  1  x  10^  flow  over  a  flat  plate  with  zero  pressure  gradient. 
This  test  case  will  be  examined  further  in  Chapters  3,  4,  and  5.  Here,  the  only  interest 
is  the  dependence  of  the  solution  on  the  modifications  to  the  turbulence  model.  For  each 
model,  the  RANS-SA  system  is  discretized  using  the  standard  weighting  scheme,  which  is 
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(a)  Original  and  modified  t]/ ji  versus  x 


(b)  Original  and  modified  Cb2X  ~  v/ versus  x 


Figure  2-2:  Comparison  of  baseline  and  modified  SA  model  diffusion  terms 
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described  in  Section  2.3.  All  results  shown  were  computed  using  p  =  3  polynomials. 

Figures  2-3  and  2-4  show  the  turbulence  model  working  variable  profiles  computed  on 
coarse  (234  elements)  and  fine  (866  elements)  meshes,  respectively.  Two  stations  are  shown: 
x/c  =  0.5  and  x/c  =  1.0.  The  profiles  are  very  similar  except  in  the  boundary  layer  edge 


(a)  x/c  =  0.5  (b)  x/c  =  0.5,  zoom 


y/c 


(c)  x/c  =  1.0 


(d)  x/c  =  1.0,  zoom 


Figure  2-3;  Turbulence  model  working  variable  profiles  for  flow  over  a  flat  plate,  computed 
using  three  versions  of  the  SA  model  on  the  coarse  mesh 

region.  In  this  region,  the  original  and  baseline  models  are  still  very  similar.  However,  while 
all  three  models  predict  D  <  0,  the  modified  model  gives  negative  f'  values  with  smaller 
magnitude.  Furthermore,  the  differences  between  the  solutions  from  the  two  models  are 
smaller  on  the  fine  mesh.  Thus,  as  the  mesh  is  refined,  the  solutions  from  the  models 
appear  to  converge. 
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(a)  xjc  =  0.5  (b)  xjc  =  0.5,  zoom 


(c)  xlc=  1.0 


(d)  x/c  =  1.0,  zoom 


Figure  2-4:  Turbulence  model  working  variable  profiles  for  flow  over  a  flat  plate,  computed 
using  three  versions  of  the  SA  model  on  the  fine  mesh 
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More  importantly,  the  differences  in  the  v  solution  have  very  little  impact  on  the  other 
flow  variables.  Figure  2-5  shows  the  velocity  profile  at  xjc  =  0.5,  computed  on  the  coarse 
mesh.  The  three  models  produce  virtually  the  same  profile.  Figure  2-6  shows  the  skin  fric- 


Figure  2-5;  Velocity  profiles  for  flow  over  a  flat  plate  at  Rex  =  5  x  10®,  computed  using 
three  versions  of  the  SA  model  on  the  coarse  mesh 

tion  distribution  for  all  three  models,  also  computed  on  the  coarse  mesh.  The  distributions 
are  practically  indistinguishable.  Furthermore,  the  drag  computed  using  the  three  models 
differs  by  only  0.02  percent.  On  the  fine  mesh,  the  spread  of  the  drag  is  less  than  5  x  10“® 
percent. 

Thus,  as  expected,  the  modifications  to  the  turbulence  model  have  little  effect  on  the 
converged  solution.  However,  as  desired,  the  negative  v  modihcations  tend  to  increase  v  in 
regions  where  i/  <  0,  leading  to  a  more  robust  scheme. 

2.3  Spatial  Discretization 

To  simplify  the  notation  for  the  discretization,  rewrite  the  RANS  equations  as 

^  +  V  •  .F(u)  -  V  •  (A(u)Vu)  =  S(u,  Vu),  (2.13) 

where  u  =  [p,  pm,  pE,  pD]'^  is  the  conservative  state  vector,  T  is  the  inviscid  flux,  AVu  is 
the  viscous  flux,  and  S  is  the  source  term. 

Let  T/j  be  a  triangulation  of  the  domain  H  C  M”"  into  non-overlapping  elements  k.  Define 
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(a)  Normal  scaling 


(b)  Log-log  scaling 


Figure  2-6;  Skin  friction  distributions  for  flow  over  a  flat  plate,  computed  using  three 
versions  of  the  SA  model  on  the  coarse  mesh 


the  function  space  by 

=  {v  G  [L\n)Y  I  V  o  /,  G  [P^KrefW,  G  Tj, 

where  r  is  the  length  of  the  state  vector,  PP  denotes  the  space  of  polynomials  of  order  p, 
and  /k  denotes  the  mapping  from  the  reference  element  to  physical  space  for  the  element 
K,  as  illustrated  in  two  dimensions  by  Figure  2-7. 


Figure  2-7:  Illustration  of  the  map  from  the  reference  element  to  the  element  k  in  physical 
space 


Then,  the  spatially  discrete  problem  is  as  follows; 


find  Uh{-,t)  G  such  that 


-L  Rh{uh,Vh) 


0,  Vv;,GV^, 


(2.14) 
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where 

-R/i(w/i,  v/i)  =  Rh,i{wh,Vh)  +  Rh,v{^h-,^h)  +  Rh,s{^h,vh), 

and  Rh,i,  Rhy-,  and  Rny  denote  the  discretizations  of  the  inviscid,  viscous,  and  source 
terms,  respectively. 


2.3.1  Inviscid  Discretization 

The  discretization  of  the  inviscid  terms  is  given  by 


where  (•)^  and  (•)“  denote  trace  values  taken  from  opposite  sides  of  a  face,  n"*"  is  the  normal 
vector  pointing  from  +  to  — ,  and  H  is  a  numerical  flux  function  for  interior  faces.  In  this 
work,  H  is  the  Roe  flux  [93].  The  inviscid  boundary  flux,  is  computed  by  evaluating 
the  flux  at  a  boundary  state,  u^,  which,  in  general,  depends  on  both  the  interior  state  and 
the  boundary  conditions.  The  boundary  condition  specification  is  discussed  in  more  detail 
in  Appendix  E. 


2.3.2  Viscous  Discretization 

The  viscous  terms  are  discretized  using  the  second  method  of  Bassi  and  Rebay  [15,  16].  To 
write  the  discretization  in  a  compact  form,  the  jump,  [-J,  and  average,  {•},  operators  are 
used.  For  interior  faces,  let 

{«}  =  +  {v}  =  + 

|s]]  =  {s'^n^  +  s~n~),  luj  =  {v^  •  •  n“). 


where  s  is  a  scalar  and  u  is  a  vector.  Then,  the  viscous  discretization  is  given  by 
Rhyi^hRh)  =  ^  /  Vv^  •  (A(w/i)Vwfe) 

KGTh 


[I 


[Whf'  ■  {A^iwh)Vvh}  -  •  ({A(w/i)V’W/i}  +  ??/{r/(w/i)})] 

(w+  -  u'’)^(A^(u')Vv+)  •  n+  +  •  n  , 
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where  denotes  the  viscous  boundary  flux  function,  Fj  and  are  auxiliary  variables, 
and  r]f  is  a  stabilization  parameter.  The  stabilization  parameter  is  set  to  r]f  =  6  for  all 
cases  shown  in  this  work.  Given  that  r]f  >  3  implies  stability  of  the  BR2  discretization  for 
triangular  meshes  (r/j  >  4  for  quadrilateral  meshes),  this  value  is  somewhat  conservative. 
The  auxiliary  variables  are  defined  by  the  following  auxiliary  problems:  for  each  face,  af, 
find  Fj  G  such  that 

Z]  /  •  {A^{wh)Th},  Wh  e  [Vir, 

for  interior  faces,  and  find  Fj  G  [V^]” 

Y,  f  rl-  = 

for  boundary  faces. 


w/  -  u 


{A^ •  n+,  yfh  G  [V^Y 


2.3.3  Source  Discretization 

Three  options  for  the  discretization  of  the  source  terms  are  considered.  The  first  option  is 
referred  to  as  the  standard  weighting  approach.  In  this  approach,  the  source  term  is  simply 
multiplied  by  the  test  function  and  integrated  over  the  domain.  Thus, 

Rh,s,sw{^h,^h)  ^  -  Y  /  VhS(wft,Vw,j).  (2.15) 

The  second  approach  is  a  modification  of  the  standard  weighting  approach  that  ensures 
the  source  term  discretization  is  dual  consistent.  The  dual  consistent  discretization  is  given 
by 

Rh,S,Dc{^h,Vh)  ^  -  Y  Vw/j) 

+  /  •  {/3(w/j,  Vw/i,  v/i)}  +  /  (u^  -  w^)/3fe(wft,  Vw/j,  Vft)  •  n. 

JTi  Jan 

The  requirements  on  the  “dual  fluxes”  /3  and  /If,  will  be  specified  in  Chapter  3,  where  a  full 
derivation  and  analysis  of  this  discretization  are  given. 

The  third  approach  is  based  on  mixed  formulations — e.g.  BR2  [15,  16]  and  LDG  [29, 
31] — where  an  auxiliary  state  variable  is  solved  for  the  gradient  of  the  state.  Specifically, 
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the  discretization  is  given  by 


-R/i,5,MF(wfe,  V/i)  =  -  ^  /  v]^S{whAh)- 


kGTh  ' 


(2.16) 


where  G  satisfies 

-  X]  /  ^h^-Th+  f  (N  •  {r/i}  +  {u}|fh]) 

+  [  uVfe-n,  Whe[Vl]"',  (2.17) 

JdQ 

and  u  =  u(w^,w^)  is  a  numerical  flux  function.  In  this  work,  u(w^,w^)  =  {w/j}  is  used. 
However,  other  valid  choices  are  available.  See  Chapter  3  for  a  complete  analysis. 

The  additional  state  variable  q/j  can  be  eliminated  by  rewriting  it  in  terms  of  Vw/j  and 
lifting  operators.  Begin  by  integrating  by  parts  on  (2.17): 


^  /  Th-Vw/j+  /  [u  -  w/jJ  •  {r/i}  +  /  {u-Wft}[[Tft]l 

Jk  JFi  JFi 

+  [  (u^-Wft)Th-n,  \fThe[Vl]''.  (2.18) 

Jan 


Define  the  lifting  operators  Vh  and  £h  (see  e.g.  [7])  by  the  following  problems;  find  rni'^h)  £ 
[V^]"^  and  ihi'^h)  £ 


-[  lu-Whj-{Th}-  [  (u^  -  w,j)Tft  •  n,  VTft  G  [V^]”  ,  (2.19) 
JTi  Jan 

-I  {u- w,j}[t/i],  VTftG[V^]''.  (2.20) 

JTi 


Then,  using  (2.18)  gives 


=  Vw/j  -  Thi^h)  -  4(wh). 


(2.21) 


Thus,  substituting  (2.21)  into  (2.16),  the  discretization  can  be  written  as 


Rh,S,MF{^h,  V/i) 


/  v^S(wft,,  Vwft,  -  r/j(w/j)  -  4(w/j)). 


2.4  Temporal  Discretization 

The  spatially  discrete  problem  (2.14)  can  be  written  as  a  system  of  ODEs.  In  particular, 
let  {vj}  for  i  =  1, . . . ,  be  a  basis  of  the  space  V^.  Then,  for  all  w/j  G  V^,  there  exists  a 
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vector  W  G  such  that 


N 

w/i(x)  =  ^  WiVi(x). 
i=l 

Thus,  the  spatially  discrete  problem  can  be  written  as  the  following  initial  value  problem; 
given  U  (0) ,  find  U  (t)  such  that 

M^  +  R{U)  =  0,  (2.22) 

at 

where  R  is  the  spatial  residual  vector, 


Ri{W)  =  Rh{wh,^i), 


and  Ai  is  the  mass  matrix. 


To  fully  discretize  the  problem,  any  number  of  standard  ODE  integration  techniques  can 
be  applied.  The  interest  of  this  work  is  steady  problems.  Thus,  the  temporal  discretization 
will  be  used  to  march  the  solution  in  time  to  a  steady  state.  In  this  case,  temporal  accuracy  is 
not  a  primary  concern.  Thus,  the  temporal  discretization  need  not  be  higher-order  accurate. 
However,  to  avoid  time  step  restrictions,  an  implicit  method  is  used.  In  particular,  the 
backward  Euler  method  is  applied.  Thus,  one  time  step  of  the  fully  discrete  problem  is  as 
follows:  given  G  ,  find  G  such  that 

-  U^)  +  R{U^+^)  =  0,  (2.23) 

where  m  is  the  iteration  number  and  At  is  the  current  time  step  length.  The  current  time 
step  is  set  based  on  an  input  CEL  number.  In  particular,  for  each  element,  define 

At,  =  CFL^, 

where  L,  is  a  measure  of  the  element  grid  spacing  and  A,  is  the  maximum  convective  wave 
speed  over  the  element.  Then, 

At  =  min  At,. 

kGTh 

At  each  successful  iteration,  the  CEL  number  is  increased  by  a  user  specified  factor,  and  a 
new  time  step  is  computed. 
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2.5  Solution  Technique 


At  each  time  step,  the  discrete,  nonlinear  problem  (2.23)  can  be  solved  using  Newton’s 
method.  For  a  solution  tolerance  of  e,  the  algorithm  is  composed  of  the  following  steps: 


1.  Set  k  =  0,V^  =  U^. 

2.  While  II {V^  -  U^)  +  R{V^)\\  >  e, 


yk+l  ^  yk 

k  =  k  +  1. 


At^  ^  dU 


3.  Set  =  V’^. 


Clearly,  this  algorithm  requires  the  solution  of  the  linear  system 


1  A.  dR 


This  system  is  inverted  using  the  GMRES  algorithm  with  element-block  ILUO  or  multigrid 
preconditioning  with  line  reordering.  See  [35,  34]  for  the  details  of  this  procedure. 

As  noted  in  Section  2.4,  steady  problems  are  the  primary  interest  of  this  work.  Thus, 
given  that  temporal  accuracy  is  not  crucial,  the  nonlinear  system  (2.23)  is  generally  not 
fully  solved  at  each  time  step.  Instead,  only  one  sub-iteration  is  taken.  In  this  case, 
is  given  by 


jjin+l  _  jjm 


1a.  dR 

Ai^^dU 


R{U^). 


Note  that  this  statement  is  not  true  when  the  unsteady  adaptation  algorithm  is  used.  For 
more  details,  see  Chapter  6. 

In  either  case,  the  solution  is  marched  forward  in  time  until  ||i?([/™)||  is  less  that  the 
desired  tolerance. 
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Chapter  3 


Dual  Consistency  Analysis 


The  focus  of  this  chapter  is  the  impact  of  dual  consistency  on  discretizations  of  source 
terms  depending  on  the  gradient  of  the  state,  like  those  seen  in  the  SA  turbulence  model 
in  Chapter  2.  Dual  consistency  provides  a  connection  between  the  continuous  and  discrete 
dual  problems.  This  concept  is  made  rigorous  in  Section  3.1,  which  gives  a  review  of  dual 
consistency  and  its  importance.  Section  3.2  analyzes  the  dual  consistency  of  the  source  term 
discretizations  shown  in  Chapter  2.  The  results  of  the  analysis  are  confirmed  by  numerical 
experiments  on  a  model  problem,  shown  in  Section  3.3.  Section  3.4  shows  results  for  the 
RANS-SA  system. 


3.1  Review  of  Dual  Consistency 

3.1.1  Definition 

Consider  the  following  primal  problem:  compute  J{u),  where  J  :  V  ^  M  is  a  functional  of 
interest,  V  is  an  appropriate  function  space,  and  tt  G  V  solves 

R{u,  v)  =  0,  Vu  G  V, 

where  i?  :  V  x  V  ^  M  is  the  weak  form  of  a  PDE  of  interest.  For  simplicity,  it  is  assumed 
that  both  the  primal  and  dual  problems  are  well-posed. 

Let  -I-  V,  where 

Vl  +  V  =  {h  =  f  +  g\fGVl,gGV}. 

Then,  consider  a  general  DC  discretization  of  the  primal  problem:  find  Uh  G  such 
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that 


Rh{uh,Vh)  =  0,  yvh£Vl, 


where  Rh  :  x  M  is  a  semi-linear  form  derived  from  the  weak  form  of  the  primal 

problem.  It  is  assumed  that  this  discrete  problem  is  also  well-posed. 

Let  Jh  ■  ^  be  the  discrete  version  of  the  functional  of  interest.  Then,  the  discrete 

dual  problem  is  given  by  the  following  statement;  find  G  such  that 

Rhl^hJivk,^)  =  Jh[uh]{vh),  Vuft  G  Vl- 

Here,  R'f^[uh]{-,'tph)  :  ^  M  is  the  linear  functional  given  by  evaluating  the  Frechet 

derivative  of  the  functional  :  Wj^  ^  M  at  Uh,  where,  for  fixed  V’ft  £ 

NxPhiwh)  =  Rhiwh.i’h),  Vrcft  G 

Similarly,  Jl,[uh]{-)  ;  ^  M  is  the  linear  functional  given  by  evaluating  the  Frechet 

derivative  of  Jh  at  Uh- 

Dual  consistency  has  been  used  in  the  analysis  of  DG  methods  by  multiple  authors.  For 
example,  [51,  50,  7]  define  dual  consistency  for  linear  problems.  A  general  definition  of  dual 
consistency  for  nonlinear  problems  is  provided  by  both  Lu  [72]  and  Hartmann  [53].  Lu  also 
defines  asymptotic  dual  consistency  for  general  nonlinear  problems.  The  definitions  used 
here  follow  [72]. 

Definition  1.  The  discretization  defined  by  the  semi-linear  form,  Rh,  and  functional,  Jh, 
is  said  to  he  dual  consistent  if,  given  exact  solutions  u  G  V  and  G  V  of  the  continuous 
primal  and  dual  problems,  respectively, 

K[u\{v,ij)  =  Jh[u]{v),  VvGWJ^. 

Definition  2.  The  discretization  defined  by  the  semi-linear  form,  Rh,  and  functional,  Jh, 
is  said  to  he  asymptotically  dual  consistent  if,  given  exact  solutions  u  G  V  and  G  V  of 
the  continuous  primal  and  dual  problems,  respectively, 

lim  I  sup  \Rh[u]{v,fi)  -  Jh[u]{v)\  ]  =  0. 

Clearly,  all  dual  consistent  discretizations  are  automatically  asymptotically  dual  consis¬ 
tent.  In  this  work,  a  discretization  will  be  referred  to  as  asymptotically  dual  consistent  only 
if  it  is  not  also  dual  consistent. 


48 


3.1.2  Importance 


For  many  types  of  discretization,  algorithms  involving  the  dual  problem  have  become  pop¬ 
ular  for  design  optimization,  error  estimation,  and  grid  adaptation  [62,  77,  36,  44,  18,  103]. 
It  is  well  known  that  dual  consistency  can  significantly  impact  the  performance  of  these 
algorithms.  For  example,  Collis  and  Heinkenschloss  [33]  showed  that  when  applying  a  dual 
inconsistent  SUPG  method  for  linear  advection-diffusion  to  an  optimal  control  problem, 
superior  results  are  obtained  using  a  direct  discretization  of  the  continuous  dual  problem 
as  opposed  to  the  discrete  dual  problem  derived  from  the  primal  discretization.  Specifi¬ 
cally,  both  the  control  function  and  the  adjoint  solution  converge  at  a  higher  rate  when  the 
continuous  dual  problem  is  discretized  directly. 

For  DG  discretizations,  Harriman  et  al.  [51,  50]  examined  symmetric  and  non-symmetric 
interior  penalty — SIPG  and  NIPG,  respectively — DG  methods  for  the  solution  of  Poisson’s 
equation.  They  showed  that  to  achieve  optimal  convergence  rates  for  a  linear  functional  out¬ 
put,  the  dual  consistent  method,  i.e.  SIPG,  must  be  used.  Lu  [72]  considered  the  impact  of 
dual  consistency  on  the  accuracy  of  functional  outputs  computed  using  DG  discretizations 
of  the  Euler  and  Navier-Stokes  equations.  He  demonstrated  the  importance  of  implementing 
the  boundary  conditions  on  the  primal  problem  in  a  dual  consistent  manner.  In  particular, 
when  using  dual  consistent  boundary  conditions,  super-convergent  functional  output  results 
are  obtained,  while,  when  using  a  dual  inconsistent  boundary  condition  treatment,  signifi¬ 
cant  degradation  of  the  output  convergence  rates  is  observed.  More  recently,  Hartmann  [53] 
has  developed  a  framework  for  analyzing  the  dual  consistency  of  DG  discretizations.  He 
uses  the  framework  to  analyze  DG  discretizations  of  many  equations,  including  the  com¬ 
pressible  Euler  and  Navier-Stokes  equations.  Eor  example,  applying  the  framework  to  the 
SIPG  discretization  of  the  Navier-Stokes  equations,  a  modification  of  the  original  SIPG 
method  to  make  the  boundary  conditions  dual  consistent  is  derived.  Similar  to  the  results 
shown  by  Lu,  Hartmann’s  modification  of  the  SIPG  scheme  produces  superior  results  to  the 
original,  dual  inconsistent  boundary  condition  treatment. 

In  addition,  it  is  well  known  that  dual  consistency  can  impact  the  convergence  of  the 
error  in  the  primal  solution  [7].  For  example,  for  many  DG  discretizations  of  Poisson’s 
equation,  standard  proofs  of  order  of  accuracy  of  the  solution  error  in  the  norm  exist. 
Typically  these  proofs  rely  on  the  Aubin-Nitsche  “duality  trick”  [96,  88]  to  obtain  an  optimal 
estimate  in  the  norm  given  an  optimal  estimate  in  the  energy  norm  [20,  7].  The  use  of  this 
duality  argument  requires  that  the  scheme  be  dual  consistent.  Thus,  some  dual  inconsistent 
methods — e.g.  NIPG  and  the  Baumann-Oden  method — do  not  achieve  optimal  accuracy 
in  the  norm,  and  dual  inconsistent  methods  that  do  achieve  optimal  accuracy  in  the 
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norm  are  typically  super-penalized  [7]. 


3.2  Analysis  of  DG  Discretizations  of  Source  Terms 

This  section  considers  DG  discretizations  of  source  terms  depending  on  the  state  and  first 
derivatives  of  the  state.  For  simplicity  of  notation,  the  analysis  will  be  done  on  a  scalar 
model  problem.  Despite  this  simplification,  the  results  easily  generalize  to  the  RANS-SA 
case. 

Let  tt  G  V  =  be  the  solution  of  the  following  scalar  problem; 

-V-{uVu)  =  f{u,Vu)  ferxeDcM*"  (3.1) 

u  =  0  for  X  G  clD, 

where  v  is  constant  and  /  G  is  the  source  term  of  interest.  As  before,  it  is  assumed 

that  this  problem  and  its  dual  are  well-posed. 

Let  J  :  V  ^  M  be  a  functional  of  interest  defined  as 

J{w)  =  [  J{w),  (3.2) 

Jn 

where  J  G  C(M).  For  simplicity,  the  functional  output  analyzed  here  does  not  include 
boundary  integrals.  The  inclusion  of  boundary  integrals  would  slightly  alter  the  boundary 
conditions  on  the  dual  problem  and  the  analysis  of  the  boundary  terms.  However,  these 
modifications  are  not  central  to  the  analysis  of  the  source  term  discretization.  For  an 
analysis  of  dual  consistency  for  equations  without  source  terms  that  includes  functionals 
with  boundary  integrals,  see  [72]. 

For  the  output  defined  in  (3.2),  the  dual  PDF  is  given  by 

—V  •  (z/VV^)  —  Dif{u,  Vu)V’  +  V  •  {D\7uf{u,  Vn)V’)  =  J'[u]  for  x  G  D,  (3.3) 

where  V’  is  the  adjoint  state,  Dif{u,Vu)  is  the  partial  derivative  of  /  with  respect  to  u 
evaluated  at  (u,  Vu)  and 

Dvufiu,  Vu)  =  [D2f{u,  Vu), . . . ,  Dn+ifiu,  Vtt)]^ 

where  Dif{u,Vu)  is  the  partial  derivative  of  /  with  respect  to  ioi  2  <  i  <  n  +  1, 

evaluated  at  (u,  Vu).  The  boundary  conditions  on  the  dual  problem  can  be  written  in  the 
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following  weak  form: 


I  dvt 


'ipq  =  0,  \/q  e  H 


The  analysis  in  Sections  3.2.1  and  3.2.2  considers  whether  the  discrete  dual  problems 
for  the  source  discretizations  defined  in  Section  2.3.3  are  consistent  with  (3.3)  and  the 
associated  boundary  conditions. 


3.2.1  The  Standard  Weighting  and  Dual  Consistent  Methods 

As  will  be  shown,  the  simple  approach  of  weighting  by  a  test  function  and  integrating  leads 
to  a  dual  inconsistent  scheme.  However,  a  dual  consistent  discretization  can  be  constructed 
by  adding  terms  to  the  discretization. 

Consider  the  following  DG  discretization:  find  Uh  G  such  that 

Rh{uh,Vh)  ^  Bh{uh,vh)  -  Vhf{uh,Vuh)  =  0,  Vu/i  G  V^,  (3.4) 

where  is  a  consistent  and  dual  consistent  bilinear  form  for  the  diffusion  operator  (e.g. 
BR2  [15]  or  LDG  [29,  31]).  Clearly,  the  source  term  discretization  in  (3.4)  is  equivalent  to 
that  given  for  the  RANS-SA  case  in  (2.15).  Furthermore,  define  the  discrete  functional  of 
interest,  Jh  '■  ^  M,  as 

Jh{wh)=  ^  /  J{wh).  (3.5) 

KdTu 

In  the  following  dual  consistency  analysis,  additional  smoothness  is  assumed  for  the  exact 
primal  and  dual  solutions,  u  and  ijj.  In  particular,  it  is  assumed  that  tt,V’  G  i/^(H).  This 
smoothness  assumption  is  common  in  the  analysis  of  dual  consistency  for  discretizations  of 
second-order  operators  [7,  72]. 

Proposition  1.  The  discretization  defined  by  the  semi-linear  form  Rh,  defined  in  (3.4), 
together  with  the  discrete  functional  Jh,  defined  in  (3.5),  is  dual  ineonsistent. 


Proof.  Linearizing  Rh  about  the  exact  solution  and  integrating  by  parts  gives 

R'h[u]iwh,Vh)  =  Bh{wh,Vh)  -  WhiDif{u,Vu)vh -V  ■  {Dyuf{u,Vu)vh)) 

KdTh 

-  I  (Nftl  •  {D^uf{u,  Vu)vh}  +  {wh}lD^uf{u,  Vu)vhj) 

JFi 


Ian 


Wh{Dyuf{u,Vu)vh)  ■  h. 


The  assumptions  V’  £  u  G  77^(H),  and  /  G  C^(M'^'''^),  imply  that  {D\!uf{u,Vu)fi’}  = 
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D\/uf{u,Vu)'ip  and  lD^uf{u,Vu)'ijj}  =  0.  Thus, 


R'h[u]{vh,'ip)  =  Bh{vh,ij)-'^  I' Vh{Dif{u,Vu)il! -V  ■  {D^uf{u,Vu)il!)) 

k&Tu 

-  [  Ivhj  ■  iDvufiu,Vu)'ilj)  -  [  Vh{DT:juf{u,Vu)'il:)-n. 

JVi  Jdn 

Evaluating  the  dual  consistency  using  the  discrete  functional  as  defined  in  (3.5)  gives 


RhbA{vh,'iP)  -  JhbAivh)  =  i^h,i{u,^)){vh)  +  {Ch,B{u,^)){vh), 


where 


i^h,i{u,i^))ivh)  ^  -  f  Ivhj- iDvuf{u,Vu)ij),  (3.6) 

JVi 

{Ch,Biu,^)){vh)  ^  -  [  VhiDyuf{u,Vu)^)-n.  (3.7) 

Jan 

In  general,  there  exists  Vh  G  such  that  at  least  {Ch,i{u,'4))){vh)  does  not  vanish.  Thus, 
the  scheme  is  dual  inconsistent.  Due  to  the  boundary  condition  on  the  dual  problem, 
the  boundary  term,  {Ch,B{ui'4^)){vh)i  will  vanish  if  {D\7uf{u,\7u)vh  ■  n)\  an  G 
However,  if  this  condition  does  not  hold,  the  boundary  term  will  also  contribute  to  the  dual 
inconsistency.  □ 


Remark  1.  It  is  possible  to  construct  a  dual  consistent  discretization  by  adding  terms  to 
the  semi-linear  form  R^.  In  particular,  define  a  new  bilinear  form, 

Rh,Dc{Wh,Vh)  =  Rh{Wh,  Vh)  +  Ahj{Wh,  Vh)  +  Ah^B{vJh,  Vh), 

where  Ahj  will  serve  to  eliminate  the  interior  face  dual  inconsistency  term,  Ci,  and  Ah^B 
will  serve  to  eliminate  the  boundary  face  dual  inconsistency,  Cb-  Furthermore,  to  maintain 
consistency,  Ahj  and  Ah^B  must  vanish  when  evaluated  at  u: 

Ah,i{u,Vh)  =  0,  \fvheVl, 

Ah,B{u,Vh)  =  0,  Vu/i  G  V^. 

The  interior  face  and  boundary  face  contributions  to  the  dual  inconsistency  are  examined 
separately.  To  eliminate  the  dual  inconsistency  from  the  interior  faces,  the  following  term 
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is  added  to  the  semi-linear  form  R^'. 


Ahj{wh,Vh)  =  /  {whi-  {i3i{wh,Vh)], 

JVi 

where  dual  consistency  requires  that  {/3j(tt,V’)}  = 

To  eliminate  the  boundary  dual  inconsistency,  the  following  term  is  added  to  Rh'- 

Ah,Biwh,Vh)  =  /  Wh0biwh,Vh)  ■  n, 

JdQ 

where  dual  consistency  requires  =  D\/uf{u,Vu)^. 


Proposition  2.  Let  be  a  dual  consistent  bilinear  form  corresponding  to  the  diffusion 
operator.  Let  be  such  that  Pi{w,  v)  =  D^uf{w,  Vw)v  for  all  w,v  ^  Let  Pb  be  such 

that  Pb{w,v)  =  D^uf{w,Vw)v  for  all  w,v  ^  Then,  the  semi-linear  form  given  by 


Rh,Dc{wh,Vh)  =  Bh{wh,Vh)  -  Vhfiwh,Vwh) 

KdTh 

+  /  Iwh}- {0i{Wh,Vh)}  +  /  Wh0b{Wh,Vh)  ■  ff 

JVi  Jan 


together  with  the  discrete  functional  Jh,  defined  in  (3.5),  is  dual  consistent. 


Proof.  Linearizing  Rh,DC  gives 


R'h,DcH{vh,'P) 


R'b[u]{vh,p)  +  [  {vhj-  {D^uf{u,Vu)p) 
JVi 


-\-  [  Vh{D^uf{u,Vu)p)  ■  n,  'ivh  G 

Jan 


Thus, 


R'h,DcbA0h,P)  -  Jh[u\{vh)  =  0,  yvh  G 


□ 


Remark  2.  The  choices  of  Pi  and  Pb  are  not  fully  determined  by  requiring  dual  consistency. 
One  valid  choice  is  given  by 


Pi{wh,Vh)  =  {Dvuf{wh,'^Wh)vh};  Pb{wh,Vh)  =  Dyuf{wh,^Wh)Vh- 
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Then, 


+ 


IWhj  ■  {Ds/uf{Wh,'^Wh)Vh}  + 


Wh{D^uf{wh,'^Wh)Vh)  ■  n. 


(3.8) 


However,  if  necessary,  one  can  construct  more  complex  functions  that  satisfy  the  dual 
consistency  requirement  as  well  as  add  stability  to  the  discretization.  For  example,  given 
that  the  dual  problem  is  a  convection-diffusion-reaction  problem,  one  option  to  add  stability 
is  to  construct  a  dual  flux  that  results  in  an  upwind  discretization  of  convective  term  of  the 
dual  problem. 


Remark  3.  In  addition  to  being  dual  consistent,  if  Bh  is  a  consistent  bilinear  form  for  the 
diffusion  operator,  the  discretization  given  in  Proposition  2  is  consistent  for  any  choice  of 
(5i  and  f3b  because,  for  the  exact  solution  u,  |u]]  =  0  and  u|ao  =  0.  Thus, 


Rh,Dc{u,Vh) 


Bh{u,Vh) 


E 


Vhfiu,Vu)  =  0, 


yvh  G  Vl 


3.2.2  The  Mixed  Formulation 


The  mixed  formulation  discretization  of  the  model  problem  is  as  follows:  find  Uh  G  such 
that 

Bh{uh,Vh)-'^  [  Vhfiuh,Vuh-rhiuh)-i'h{uh))  =  0,  G  V^,  (3.9) 

where  fh  and  ih  are  lifting  operators  analogous  to  those  defined  in  Section  2.3.3.  In  par¬ 
ticular,  they  are  defined  by  the  following  problems;  find  fh{uh)  G  [V^]”  and  G 

such  that 


Uhj  •  {Th} 


-  /  {n{uh)  -  Uh)Th  ■  n, 
JdQ 

-  /  {u{Uh)  -  UhWhj, 

JTi 


Vr.G[V,-]^ 

Vr.G[V,^]r 


(3.10) 

(3.11) 


To  complete  the  scheme,  one  must  define  the  numerical  flux  functions  u  and  For  the 
remainder  of  the  chapter,  it  is  assumed  that  these  fluxes  have  the  following  properties; 
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1.  for  each  interior  edge,  e,  there  exists  a  constant  vector  de  such  that 

u{wh)  =  {Wh}  +  de  ■ 


2.  =  0. 


Remark  4.  While  the  assumptions  on  the  numerical  fluxes  may  seem  restrictive,  to  the  best 
of  the  author’s  knowledge,  all  existing  mixed  formulations  that  are  consistent  use  fluxes  sat¬ 
isfying  these  assumptions  for  problems  with  homogeneous  Dirichlet  boundary  conditions  [7]. 

The  dual  consistency  of  this  discretization  is  considered  in  two  parts.  First,  a  restrictive 
condition  is  assumed  which  implies  that  the  scheme  is  dual  consistent.  Then,  the  proof  is 
extended  to  the  general  case,  where  only  asymptotic  dual  consistency  can  be  shown. 


Proposition  3.  Let  D\/uf{u,Vu)‘ip  G  Then,  the  DG  formulation  defined  in  (3.9) 

together  with  the  functional  Jh  defined  in  (3.5)  forms  a  dual  consistent  discretization. 


Proof.  Noting  that  the  lifting  operators  and  ih  are  linear  functionals  and  that  rh{u) 
£h{u)  =  0,  linearizing  Rh  about  the  exact  solution  gives 

Rh[u]{wh,Vh)  =  Bh{wh,Vh)-'^  /  VhDif{u,Vu)wh 

KGTh 

/  VhD^uf{u,Vu)-{Vwh-rh{wh)-ih{vJh))- 

Thus,  integrating  by  parts  gives 

R'ffu\{wh,Vh)  =  Bh{wh,Vh)-'^  Wh{Dif{u,Vu)vh-V-{Dx7uf{u,Vu)vh)) 

-  [  (Iw'J  •  {Dvuf{u,  Vu)vh}  +  {wh}lD^uf{u,  Vu)i;fel) 


idn 


Wh{B>\7uf{u,Vu)vh)  ■  h 


+  ^  VhDvufiu,Vu)  ■  (fhiwh)  +ihiwh))- 


(3.12) 


Using  the  assumptions  on  u  and  combined  with  the  hypothesis  that  Ds/uf{u,  G 
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Kl", 

^  f  iD’vufiu,'^u)ip)-fh(wh)  ^  /  Iwft])  ■  {Dv„/(u,  Vu)?/!} 

KGTh 

+  /  Wh{D^ufiu,Vu)ip)-n, 

Jdn 

f  {D^uf{u,Vu)'il))  ■  ih{wh)  =  -  f  {de-lwh})lD^uf{u,Vu)i;}. 

K&Th  ^ 

Substituting  (3.13)  and  (3.14)  into  (3.12)  gives 


(3.13) 

(3.14) 


R'h[u\{wh,Vh)  =  Bh{wh,Vh)-'^  Wh{Dif{u,Vu)vh-V-{Dx7uf{u,Vu)vh)) 

„  ^  rp  J  K, 


KGTh  ' 

{u{wh)lDs7uf{u,Vu)vh}) . 


Furthermore,  by  the  assumptions  on  the  smoothness  of  u,  V’,  and  /,  lD\/uf{u,Vu)^}  =  0. 
Thus, 

R'hW]{vh,'ip)  =  Bhivh,ij)  -  ^  f  VhiDif{u,Vu)ij  -  V  ■  {D^uf{u,Vu)il!)) . 

KGTh 

Finally, 

R^h[u]ivh,^P)  -  Ji[u]{vh)  =  0,  yvh  e  Wl, 

which  implies  that  the  scheme  is  dual  consistent.  □ 

Of  course,  in  general,  the  assumption  of  Tlvn/('a,  Vtt)'0  G  is  not  realistic.  For 

Dvuf{u,Vu)^  0  Proposition  3  does  not  hold.  In  this  case,  the  mixed  formulation  is 

only  asymptotically  dual  consistent. 

Proposition  4.  If  D\ruf{u,Vu)'ip  G  {!})]"• ,  where  I  <  k  <  p,  then  the  DG  diseretiza- 

tion  defined  in  (3.9)  together  with  the  functional  Jh  defined  in  (3.5)  forms  an  asymptotically 
dual  consistent  discretization. 

To  simplify  the  proof,  some  preliminary  notation  and  lemmas  are  required.  In  particular, 
let  £h  denote  the  set  of  all  faces  in  the  triangulation  Th-  Define  the  jump  operator,  !•],  on 
boundary  faces  by  |s]]  =  sn  for  scalar  quantities  and  \v\=  v-n  for  vector  quantities.  Define 
the  average  operator,  {•},  on  boundary  faces  by  {u}  =  v. 

For  the  following  lemmas,  it  is  assumed  that  =  (V^  +  n  Hq{Q).  Further¬ 

more,  it  is  assumed  that  the  set  of  triangulations,  [T^lh^Q,  is  quasi-uniform  (see  [38,  88]  for 
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definition). 

Lemma  1.  There  exists  a  norm,  |||  •  |||*  :  ^  M,  and  a  constant,  c,  such  that 


h  IINIIo,e  <  Y  ^  clllvlll*,  'iv  G  Wl, 

e&Sh  e&Sh 

where  Ke  is  such  that  e  C  due,  h  =  max^eTh  h^,  and  h^  =  \x  —  y\. 

Proof.  An  example  of  such  a  norm  is  used  by  Arnold  et  al.  [7].  In  particular, 

llblll*  =  Y  +  hl\v\Y)  +  Y  lke(M)llo,0> 

K^Th 


where 

/r.(W).f=- /'hUt}  Vfelva” 

JQ.  J  e 

For  the  proof  of  the  lemma  for  this  norm,  see  [7],  Section  4.1,  or  [20],  Lemma  2.  Note  that 
only  the  existence  of  such  a  norm,  not  its  particular  form,  is  important  here.  □ 

Lemma  2.  For  a  face,  e  G  Sh,  such  that  e  C  dn,  there  exists  a  constant,  c,  such  that,  for 
all  V  G  H^{k)  and  w  G  L^(e), 

vw\  <  c/i“^/^(||u||o,K  +  h«,|u|i,«,)||u;||o,e- 


Proof.  Apply  the  Cauchy-Schwarz  inequality,  and  then  use 

||^'||o,e  <  c/i“^/^(||u||o,«;  +  Vu  G 

which  is  a  standard  trace  theorem  [7].  □ 

Lemma  3.  For  all  w  G  where  1  <  k  <  p,  there  exists  a  constant,  c,  such  that 

^  ||u;  -  n)](u;)||?  <  ch^\w\k+i,Q, 

\KGTh  / 

where  II)]  :  L^(n)  ^  is  the  L?' {Ft) -orthogonal  projection  onto  V^. 

Proof.  If  IIk  :  L^(k)  ^  P^  is  the  L^(K)-orthogonal  projection  onto  P^,  then 

nP(u)U  =  np(uU),  yvGL\n). 
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To  complete  the  proof,  apply  Proposition  1.134(iii)  from  [38]  to  each  element  n  and  sum 
over  the  elements.  □ 


Proof  (Proposition  4)-  Define  tt  G  [V^]”  by 

=  nft((Dv«/(u,  Vu)'4))j),  for  j  =  1, . . . ,  n. 

Furthermore,  define  e  G  [L^(D)]"'  by 

ej  =  {D^uf{u,Vu)'4!)j  -  TTj,  for  j  =  1, . . . ,  n. 

By  assumption,  D\7uf{u,Vu)il)  G  Thus,  ej^  G  [H^{k)]^,  for  all  k  G  T^.  From  the 

proof  of  Proposition  3,  it  is  clear  that 


Ehivh)  ^  R'hiuKvh,^^  -  Jh[u]{vh)  =  ^  f  iDyufiu,Vu)f;)  ■  (Ph{vh) - 

K&Th 

-  Vu)V'),  yvh  G 

where  F  =  Fj  U  dQ.  Thus, 


Eh{vh)  =  [{^  +  ^ 

K&Th 

-  jlvhl-{T^  +  e),  Mvh^yVl. 
From  (2.19),  (2.20),  and  the  assumptions  on  u  and  vf, 

^  [^-rhivh)  =  [lvhj-{Tr}, 

^  f  ^•4(^'h)  =  -  f  {de  ■  IVhM^j- 

cC.'T,  ^ 


(3.15) 


(3.16) 

(3.17) 


K&Th 


Substituting  (3.16)  and  (3.17)  into  (3.15)  gives 

Eh{vh)  =  '^  [  d-  {rh{vh)  +  h{vh))  -  [  bh}  •  -  [  (de  • 

By  the  definition  of  e, 
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Furthermore,  since  lD\/uf{u,Vu)^}  =  0,  it  is  clear  that  Ivf]]  =  — Thus, 


Eh{vh)  =  -  [  Ivhj  •  {e}  +  [  {de  ■  NDH- 
Jr  JTi 

Then,  using  the  assumption  on  u,  Eh{vh)  can  be  rewritten  as 


E, 


.(vh)  =  -  Vhe-  n-  lufej 
Jdn  JTi 


^  -  cTe  •  n+'j  +  Q  +  cTe  •  )  r 


where  (•)'*'  and  (•)  refer  to  trace  values  taken  from  opposite  sides  of  an  interior  face,  and 
n"*"  is  the  outward  pointing  unit  normal  for  k~^. 

Applying  Lemma  2  to  each  edge  in  the  triangulation,  one  can  show  that 


\Eh{Vh)\  ||Klllo,e 

eeSh i=i 


Thus, 


n 

mn)\  <  EEh''-T{ikjiiL.)‘''"  iihiiio. 


< 


eeSh J=1 

n 

EE 

eG£/i  i=l 


1/2 


K&Th 


Milo, 


n 

■  ,  ,  1/2- 

- 

s  GE 

j  E  ii'jiiL) 

X 

E II 

i=i 

_ee£h 

Finally,  applying  Lemmas  1  and  3  gives 


\Eh{vh)\  <  C\\\vh\\\*h^'^\{D^uf{u,Vu)'4))j\k+i^Q_. 

/=i 

Thus,  as  h  ^  0,  \Eh{vh)\  0  for  all  Vh  G  VV^,  which  implies  that  the  scheme  is  asymptoti¬ 
cally  dual  consistent.  □ 


3.3  Model  Problem  Results 

As  a  demonstration  of  the  effects  of  dual  consistency,  a  simple  test  problem  based  on  a 
nonlinear  ODE  is  considered.  The  effect  of  dual  consistency  on  the  convergence  rates  of  the 
solution  and  adjoint  solution  errors  as  well  as  a  simple  functional  output  is  demonstrated. 
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Consider  the  following  ODE: 


-{{ly  +  u)ua;)x  -  cUxUx  =  g  for  XG  (0,1), 
u(0)  =  u{l)  =  0, 

where  u  =  1  and  c  =  ^ .  Setting 

g{x)  =  7r^((i^  +  sin(7rx))  sin(7rx)  —  (1  +  c)  cos^(7rx)), 
it  is  easy  to  show  that  the  exact  solution  is  given  by 

u(x)  =  sin(7rx). 

This  nonlinear  problem  has  been  solved  using  three  discretizations:  the  standard  weight¬ 
ing  method  as  shown  in  (3.4),  a  dual  consistent  method  with  a  penalty  term  as  shown  in 
(3.8),  and  an  asymptotically  dual  consistent  mixed  method  with  u  =  {tt}.  In  all  cases,  the 
BR2  scheme  is  used  to  discretize  the  nonlinear  diffusion  operator. 

Figure  3-1  shows  the  error  in  the  primal  solution  versus  grid  refinement.  The  error  is 
measured  in  a  broken  norm  defined  by 

KdTh 

In  this  norm,  all  three  schemes  produce  essentially  the  same  error  in  the  primal  solution. 
However,  as  shown  in  Figure  3-2,  the  dual  consistent  and  asymptotically  dual  consistent 
schemes  produce  superior  results  when  the  error  is  measured  in  the  norm.  In  particular, 
for  the  dual  inconsistent  discretization,  the  norm  of  the  error  is  proportional  to  hP 
for  even  p  and  proportional  to  for  odd  p.  It  is  not  clear  why  the  even  and  odd  p 
results  show  different  asymptotic  rates,  but  similar  results  have  been  observed  for  other 
dual  inconsistent  discretizations  [51]. 

Alternatively,  the  dual  consistent  and  asymptotically  dual  consistent  discretizations  give 
0{hP~^^)  for  all  p  tested.  Furthermore,  it  is  interesting  to  note  that  the  asymptotically 
dual  consistent  method  produces  essentially  exactly  the  same  results  as  the  dual  consistent 
discretization.  This  fact  shows  that,  for  the  asymptotically  dual  consistent  scheme,  the 
contribution  of  the  dual  inconsistency  error  to  the  total  error  is  sufficiently  high-order  that 
it  does  not  degrade  the  asymptotic  rate  convergence  rate  of  the  error.  This  result  is  not 
surprising  given  the  form  of  the  dual  consistency  error  derived  in  Section  3.2.2. 
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lOQjfhg/h) 


(a)  Dual  consistent 


(b)  Standard  weighting 


(c)  Mixed  formulation 


Figure  3-1;  Primal  error  in  the  broken  norm  for  the  scalar  model  problem 
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L  (!2) 


(a)  Dual  consistent 


(b)  Standard  weighting 


(c)  Mixed  formulation 


Figure  3-2:  Primal  error  in  the  norm  for  the  scalar  model  problem 
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Let  the  output  of  interest  be 


where  w{x)  =  2sin(7rx).  Then,  examining  the  adjoint  solution  error,  one  can  see  that  the 
dual  consistent  and  asymptotically  dual  consistent  schemes  are  superior  for  computing  the 
adjoint.  Figure  3-3  shows  the  adjoint  error  in  the  broken  norm.  The  adjoint  error 


^  -  r1  - 4 - < 

-e-p=2 

-<-p=3 

-^p=4 

> 

(a)  Dual  consistent 


(b)  Standard  weighting 


(c)  Mixed  formulation 


Figure  3-3:  Adjoint  error  in  the  broken  norm  for  the  scalar  model  problem 

is  computed  relative  to  a  40th  order  solution  of  a  Galerkin  spectral  discretization  of  the 
dual  problem.  When  using  the  dnal  inconsistent  discretization,  the  broken  norm  of  the 
adjoint  error  does  not  converge  to  zero  with  grid  refinement.  For  the  dual  consistent  and 
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asymptotically  dual  consistent  schemes,  this  error  converges  at  0{h^).  Similarly,  Figure  3- 
4  shows  that  the  norm  of  the  adjoint  error  converges  at  0{h)  when  using  the  dual 
inconsistent  scheme,  regardless  of  p,  while,  for  the  dual  consistent  and  asymptotically  dual 
consistent  schemes,  this  error  converges  at 
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(a)  Dual  consistent 


(b)  Standard  weighting 


(c)  Mixed  formulation 


Figure  3-4;  Adjoint  error  in  the  L?'  norm  for  the  scalar  model  problem 

Since  the  exact  solution  is  known,  computing  the  exact  functional  output  is  trivial, 
enabling  comparison  of  the  computed  result  with  the  exact  value,  J{u)  =  1/4.  Figure  3-5 
shows  the  error  in  the  computed  functional  for  the  three  discretizations  considered.  The 
figure  shows  that,  as  in  the  state  and  adjoint  results,  the  performance  of  the  dual  consistent 
and  asymptotically  dual  consistent  schemes  is  very  similar.  Both  schemes  achieve  for 
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|J(U)  -  J(u.  )| 


(a)  Dual  consistent 


(b)  Standard  weighting 


(c)  Mixed  formulation 


Figure  3-5:  Functional  output  error  for  the  scalar  model  problem 
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1  <  p  <  4.  However,  for  the  dual  inconsistent  scheme,  the  convergence  rate  of  the  functional 
is  0{hP)  for  even  p  and  for  odd  p.  Thus,  the  dual  consistent  and  asymptotically 

dual  consistent  discretizations  predict  the  functional  with  greater  accuracy  than  the  dual 
inconsistent  discretization  for  similar  numbers  of  degrees  of  freedom. 

3.4  RANS  Results 

Three  test  cases  are  considered.  First,  the  adjoint  solution  and  its  derivatives  are  examined 
for  flow  over  a  flat  plate  solved  using  all  three  source  term  discretizations.  Second,  a  grid 
refinement  study  is  conducted  for  flow  around  a  NACA  0012  airfoil  using  both  the  standard 
weighting  and  dual  consistent  discretizations.  Finally,  a  RAF  2822  airfoil  case  is  solved 
using  the  dual  consistent  scheme. 

3.4.1  Flat  Plate  Adjoint  Results 

To  begin,  adjoint  solutions  for  flow  over  a  flat  plate  are  examined.  In  particular,  the  test 
case  is  M^o  =  0.25,  Rec  =  1  x  10^  flow  over  a  flat  plate  with  zero  pressure  gradient.  The 
output  of  interest  is  the  drag  force  on  the  plate. 

Figures  3-6  and  3-7  show  p  =  3  adjoint  solution  profiles  computed  on  an  866  element 
mesh  and  on  a  3464  element  mesh  generated  by  uniformly  refining  the  866  element  mesh. 
Even  on  the  coarser  mesh,  it  is  difficult  to  see  a  difference  between  the  profiles.  This  fact 
is  consistent  with  the  adjoint  results  from  the  model  problem.  Specifically,  for  the  model 
problem,  the  adjoint  error  measured  in  the  norm  converges  for  all  three  schemes  but  at 
different  rates. 

However,  as  can  be  seen  in  Figures  3-8  and  3-9,  there  are  clearly  differences  in  the 
derivatives  of  the  adjoint.  In  particular,  the  figures  show  that  the  y  derivatives  of  the 
adjoint  computed  using  the  standard  weighting  scheme  are  highly  oscillatory,  while  the 
dual  consistent  and  mixed  formulation  results  are  relatively  smooth.  Furthermore,  the 
oscillations  in  the  derivatives  produced  by  the  standard  weighting  scheme  are  not  decreasing 
with  grid  refinement.  This  result  is  consistent  with  the  results  obtained  for  the  adjoint  error 
measured  in  the  broken  norm  for  the  model  problem.  In  that  case,  the  adjoint  error  did 
not  converge  for  the  standard  weighting  scheme.  Finally,  the  figure  shows  that,  while  the 
dual  inconsistency  in  the  standard  weighting  scheme  is  caused  only  by  the  discretization  of 
the  turbulence  model — none  of  the  other  equations  have  gradient  dependent  source  terms — 
it  adversely  affects  the  adjoint  components  corresponding  to  the  other  equations  as  well. 
This  adverse  effect  is  not  surprising  given  the  coupling  between  the  turbulence  model  and 
the  flow  equations. 
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y/c 


(a)  866  element  mesh 


y/c 


(b)  3464  element  mesh 


Figure  3-6:  Comparison  of  x-momentum  adjoint  solution  profiles  for  dual  consistent,  mixed 
formulation,  and  standard  weighting  discretizations  at  x/c  =  0.5  for  flow  over  a  flat  plate 
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(a)  866  element  mesh 


(b)  3464  element  mesh 


Figure  3-7:  Comparison  of  turbulence  model  adjoint  solution  profiles  for  dual  consistent, 
mixed  formulation,  and  standard  weighting  discretizations  at  x/c  =  0.5 


68 


(a)  866  element  mesh 


(b)  3464  element  mesh 


Figure  3-8:  Comparison  of  x-momentum  adjoint  y-derivative  profiles  for  dual  consistent, 
mixed  formulation,  and  standard  weighting  discretizations  at  x/c  =  0.5  for  flow  over  a  flat 
plate 
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y/c 


(a)  866  element  mesh 


y/c 


(b)  3464  element  mesh 


Figure  3-9:  Comparison  of  turbulence  model  adjoint  y-derivative  profiles  for  dual  consistent, 
mixed  formulation,  and  standard  weighting  discretizations  at  x/c  =  0.5  for  flow  over  a  flat 
plate 
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3.4.2  NACA  0012  Refinement  Stndy 


To  examine  the  accuracy  of  the  discretizations,  a  grid  refinement  study  is  conducted  for 
-^oo  =  0.25,  Rcc  =  1  X  10^,  a  =  0  flow  over  a  NACA  0012  airfoil.  The  airfoil  is  modified 
slightly  such  that  the  trailing  edge  has  zero  thickness.  Specifically,  the  airfoil  is  defined  by 

y  =  ±0.6(0.2969Vx  -  0.1260x  -  O.SSlOx^  +  0.2843x^  -  O.lOSOx"^), 

for  0  <  X  <  1.  The  outer  boundary  of  the  computational  domain  is  the  circle  of  radius 
40,  centered  at  the  origin.  To  generate  a  family  of  meshes,  a  coarse,  1280  element  linear 
mesh  of  this  domain  was  generated  by  meshing  a  rectangle  and  then  mapping  the  mesh 
nodes  to  the  physical  domain  using  transfinite  interpolation.  The  resulting  mesh  is  shown 
in  Figure  3-10.  Given  the  coarse  mesh,  medium,  fine,  and  ultra-fine  linear  meshes  were 
generated  by  uniformly  refining  the  coarse  mesh. 

Results  are  presented  for  p  =  1,2,  3, 4  for  the  dual  consistent  and  standard  weighting 
discretizations  on  the  coarse,  medium,  and  fine  meshes.  Only  p  =  1,2  results  have  been 
obtained  on  the  ultra-fine  mesh,  and  no  mixed  formulation  results  are  presented  because 
this  discretization  has  not  been  implemented  in  parallel. 

All  results  are  isoparametric.  Thus,  curved  meshes  are  required.  To  generate  these 
curved  meshes,  a  linear  elasticity  based  node  movement  algorithm,  described  in  Section  5.1.3, 
was  applied  to  the  linear  meshes.  To  illustrate  the  differences  between  the  linear  and  curved 
meshes,  the  q  =  1  and  q  =  3 — where  q  denotes  the  polynomial  order  of  the  mesh — versions 
of  the  coarse  mesh  are  shown  in  Figure  3-11.  The  q  =  3  mesh  is  plotted  by  linearly  inter¬ 
polating  between  the  edge  nodes.  The  figure  shows  the  leading  edge  region  of  the  airfoil, 
where  the  differences  between  the  q  =  1  and  q  =  3  meshes  are  the  largest.  Far  away  from 
the  airfoil,  the  two  meshes  are  essentially  the  same. 

To  examine  the  drag  error,  an  exact  value  of  the  drag  must  be  chosen.  Since  the  exact 
solution  is  not  known,  the  “exact”  drag  in  this  work  is  computed  by  extrapolating  high- 
order  results.  Specifically,  the  error  is  computed  relative  to  a  drag  value  of  75.612866  counts, 
which  is  obtained  by  extrapolating  the  dual  consistent,  p  =  5  drag  results  from  the  coarse 
and  medium  meshes,  assuming  the  drag  error  for  p  =  5  is  proportional  to  h^.  This  rate 
was  chosen  after  examining  the  order  of  accuracy  obtained  by  the  p  =  1,2,  3, 4  results  for 
different  p  =  5  rate  assumptions.  The  results  of  this  examination  are  shown  in  Figure  3-12. 
The  figure  shows  the  order  of  accuracy  obtained  for  p  =  1, 2,  3, 4  by  refining  from  the  coarse 
mesh  to  the  fine  mesh  versus  the  assumed  order  of  accuracy  used  to  extrapolate  the  p  =  5 
drag  results  to  the  “exact”  value.  The  p  =  1,2  results  are  relatively  insensitive  to  the 
p  =  5  rate  assumption.  This  result  gives  confidence  that  the  computed  error  and  order  of 
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(a)  Full  domain 


(b)  Zoom 


Figure  3-10:  Coarse  (1280  elements),  linear  mesh  of  the  NACA  0012  airfoil 


(a)  g  =  1 


(b)  g  =  3 


Figure  3-11:  Leading  edge  of  the  q  =  1  and  q  =  3  versions  of  the  coarse  mesh  of  the  NACA 
0012  airfoil 


Figure  3-12:  Order  of  accuracy  obtained  for  p  =  1,  2, 3, 4  versus  order  of  accuracy  assumed 
for  p  =  5  for  drag  error  for  flow  over  a  NACA  0012  airfoil,  computed  using  the  dual 
consistent  discretization 

accuracy  for  p  =  1,2  are  very  close  to  the  true  values.  Alternatively,  the  p  =  3,4  results 
are  quite  sensitive  to  the  assumed  rate.  However,  it  appears  that  assumed  rates  between 
two  and  three  give  the  most  consistent  results.  Assuming  p  =  5  is  converging  faster  than 
approximately  0{h^)  leads  to  a  scenario  where  p  =  3  obtains  between  0{h^)  and  0{h‘^) 
while  p  =  4  gives  between  0{h?‘)  and  0{h^).  This  result  makes  little  sense.  On  the  other 
hand,  assuming  p  =  5  obtains  0{h^)  gives  that  p  =  3,4  both  achieve  roughly  0{h^).  While 
this  rate  is  suboptimal,  is  can  be  reasonably  explained  if  the  solution  is  not  smooth.  Thus, 
0{h^)  is  chosen  here. 

Figure  3-13  shows  the  drag  error  results  obtained  using  the  dual  consistent  discretiza¬ 
tion.  In  particular.  Figure  3-13(a)  shows  the  drag  error  versus  grid  refinement.  The  p  =  1 
order  of  accuracy  is  slightly  better  than  the  expected  optimal  rate  of  0{h?‘),  and  the  p  =  2 
accuracy  appears  optimal  at  approximately  0{h‘^).  Alternatively,  p  =  3,4  do  not  achieve 
the  asymptotic  convergence  rates  expected  for  smooth  problems.  While  the  precise  rates 
obtained  are  dependent  on  the  assumed  p  =  5  rate,  no  assumed  rate  gives  optimal  perfor¬ 
mance  for  both  p  =  3  and  p  =  4,  as  shown  in  Figure  3-12.  This  suboptimal  behavior  is 
attributed  to  under-resolution  in  the  boundary  layer  edge  region.  Neglecting  the  effect  of 
the  viscosity  of  the  fluid,  the  RANS-SA  system  has  discontinuous  first  derivatives  at  the 
boundary  layer  edge  in  both  the  primal  and  adjoint  solutions.  The  discontinuities  in  the 
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Iog2(ho/h) 


(a)  Drag  error  versus  h  for  p  —  1  through  p  =  5  (p  =  1  :  p  =  2  :  o, 

p  =  3:*,p  =  4:<l,  p  =  5:  x) 


(b)  Drag  error  versus  DOF 


Figure  3-13:  Drag  error  for  flow  over  a  NACA  0012  computed  using  the  dual  consistent 
discretization 
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adjoint  derivatives  can  be  seen  in  Figures  3-8  and  3-9.  While  the  addition  of  viscosity  for¬ 
mally  eliminates  these  discontinuities,  given  that  the  solution  in  this  region  is  not  resolved, 
the  effect  on  the  numerical  solution  is  the  same. 

Figure  3-13(b)  shows  the  drag  error  results  again.  In  this  figure,  the  drag  error  versus 
DOF  is  shown  for  p  refinement  on  the  coarse,  medium,  fine,  and  ultra-fine  meshes.  The  figure 
shows  that,  on  each  mesh,  as  p  is  increased,  the  drag  error  eventually  becomes  proportional 
to  DOF~^/^ .  In  this  regime,  the  error  is  dominated  by  the  boundary  layer  edge  region,  as 
will  be  confirmed  by  examining  an  error  estimate  in  Section  4.3.2.  To  further  confirm  that 
the  boundary  layer  edge  singularity  is  the  root  cause  of  the  suboptimal  accuracy,  a  laminar 
case  is  examined  using  the  same  set  of  meshes  used  here.  The  laminar  case  has  no  boundary 
layer  edge  singularity,  and,  as  expected,  the  order  of  accuracy  does  not  asymptote  to  0{h^) 
as  the  resolution  increases.  The  results  of  this  study  are  shown  in  Appendix  B. 

The  results  of  the  refinement  study  for  the  standard  weighting  discretization  are  shown 
in  Figure  3-14.  As  in  the  dual  consistent  case,  the  p  =  1  discretization  gives  better  than 
expected  results.  However,  the  effect  of  the  dual  inconsistency  can  be  observed  for  the  p  =  2 
scheme,  which  obtains  only  0{h?)  accuracy.  This  rate  is  consistent  with  the  suboptimal 
results  observed  for  the  standard  weighting  discretization  applied  to  the  model  problem. 
For  p  >  3,  it  appears  that  errors  due  to  the  boundary  layer  edge  singularity  are  dominant. 
Thus,  no  degradation  of  the  accuracy  due  to  the  dual  inconsistency  is  observed.  In  fact,  in 
this  regime,  the  standard  weighting  discretization  appears  to  be  slightly  more  accurate  for 
this  case. 

3.4.3  RAE  2822  Refinement  Study 

To  demonstrate  the  performance  of  the  dual  consistent  scheme  on  a  slightly  more  difficult 
problem,  a  RAE  2822  airfoil  case  is  considered.  The  freestream  flow  conditions  are  M^o  = 
0.6,  Rcc  =  6.3  X  10®,  a  =  2.57°.  The  computational  domain  is  the  circle  of  radius  40, 
centered  at  the  origin,  which  is  the  leading  edge  of  the  airfoil.  The  meshes  are  generated 
analogously  to  those  used  for  the  NACA  0012  case.  That  is,  a  coarse,  linear  mesh  containing 
4000  elements  was  generated.  Figure  3-15  shows  this  coarse  mesh.  The  linear  coarse  mesh 
was  uniformly  refined  twice  to  produce  meshes  of  16000  and  64000  elements.  Then,  high- 
order  meshes  were  constructed  using  the  linear  elasticity  node  movement  scheme.  Results 
are  presented  for  p  =  1, 2,  3, 4,  5  for  the  coarse  and  medium  meshes.  For  the  fine  mesh,  only 
p  =  1,2,3  solutions  have  been  computed. 

Figure  3-16  shows  the  drag  error  results  obtained  using  the  dual  consistent  discretization. 
The  error  is  computed  relative  to  an  “exact”  drag  value,  which  is  computed  by  extrapolating 
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Iog2(ho/h) 


(a)  Drag  error  versus  h  for  p  —  1  through  p  =  5  (p  =  1  :  p  =  2  :  o, 

p  =  3:*,p  =  4:<l,  p  =  5:  x) 


(b)  Drag  error  versus  DOF 


Figure  3-14:  Drag  error  for  flow  over  a  NACA  0012  computed  using  the  standard  weighting 
discretization 
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the  p  =  5  drag  results,  assuming  that  the  error  is  proportional  to  h^.  This  rate  is  chosen  to 
be  consistent  with  the  NACA  0012  case.  The  resulting  value  is  94.542604  counts. 

The  figure  shows  the  drag  error  versus  grid  refinement  and  DOF.  The  p  =  1  order  of 
accuracy  is  better  than  expected  at  approximately  0{h^),  and  the  p  =  2  accuracy  oscillates 
slightly  about  the  optimal  0(/i^).  However,  as  in  the  NACA  0012  case,  the  p  =  3,4  results 
do  not  achieve  optimal  accuracy.  Unlike  the  NACA  0012  case,  the  suboptimal  accuracy 
here  is  attributed  to  under-resolution  of  the  trailing  edge  region.  This  under-resolution  will 
be  demonstrated  by  examining  elemental  error  estimates,  as  shown  in  Section  4.3.3. 
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(a)  Full  domain 


(b)  Zoom 


Figure  3-15;  Coarse  (4000  elements),  linear  mesh  of  the  RAE  2822  airfoil 
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(a)  Drag  error  versus  h  for  p  =  1  through  p  =  5  (p  =  1  :  p  =  2  :  o, 

p  =  3:*,p  =  4:<l,  p  =  5:  x) 


(b)  Drag  error  versus  DOF 


Figure  3-16:  Drag  error  for  flow  over  a  RAE  2822  computed  using  the  dual  consistent 
discretization 
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Chapter  4 


Output  Error  Estimation 


This  chapter  details  the  error  estimation  approach  applied  in  this  work.  Sections  4.1  and  4.2 
show  the  error  estimation  algorithm,  which  is  based  on  the  method  of  Dual  Weighted 
Residuals,  due  to  Becker  and  Rannacher  [17,  18],  and  draws  heavily  on  recent  work  by 
Fidkowski  and  Darmofal  [42,  39].  The  main  contribution  of  this  chapter  is  the  application 
of  the  algorithm  to  high-order  discretizations  of  the  RANS  equations  on  curved,  boundary 
conforming  meshes.  Numerical  results,  shown  in  Section  4.3,  demonstrate  that  the  error 
estimation  algorithm  produces  acceptable  estimates  of  the  drag  error  for  three  high  Re 
problems. 

4.1  Motivation 

Let  Rh{-,  •)  be  a  semi-linear  form,  and  let  uh  £  Vh  be  such  that 

Rh(u_h-,  v/f)  =  0,  VvH  G  Vh,  (4.1) 

where  Vh  is  an  appropriate  finite  dimensional  function  space — e.g.  V^  as  defined  in  Sec¬ 
tion  2.3.  Furthermore,  let  u  G  V  be  the  exact  solution  of  the  PDF  of  interest.  Then,  for  a 
general  output  of  interest,  J{  ),  define  the  following  dual  problem:  find  '0  G  V  such  that 

.Rh(u,uh;  V,  V’)  =  77(u,uh;  v),  VvgVh  +  V,  (4.2) 

where  Rh  and  J  are  mean-value  linearizations  given  by 

.Rh(u,  uh;  v,’w)  =  f  Rff[eu+ {1  -  e)uH]i^,w)de, 

Jo 

77(u,uh;v)  =  [  J'[0u+ {1  -  0)uH]{v)d6. 

Jo 
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Then, 


RHiu,UH;u- UH,w)  =  Rh{u,w)  -  Rh{uh,w), 

J{u,uh;u-uh)  =  J{u)-J{uh). 

Assuming  that  Rh{u,-w)  =  0,  Vw  G  Vh  +  V  and  using  (4.1)  and  (4.2),  the  error  in  the 
output  may  be  written  as 

J^(u)  -  J{uh)  =  -Rh{uh,  ^  (4-3) 

for  arbitrary  G  Vh-  Alternatively,  by  defining  the  adjoint  residual, 

.R|)(u,  Uh;  v,w)  =  RH{u,UH]^r,w)  -  J(u,u/^-;  v), 

it  is  possible  to  express  the  error  in  terms  of  the  exact  primal  solution.  In  particular,  one 
can  show  that 

J^(u)  -  J{uh)  =  -Rh{u,  uh;  u  -  UH,  (4-4) 

again  for  any  ipH  ^ 


4.2  Implementation 

Given  that  (4.3)  and  (4.4)  depend  on  the  unknown  exact  solutions,  it  is  necessary  to  make 
some  approximations  to  compute  error  estimates.  To  begin,  all  mean-value  linearizations 
are  replaced  by  linearizations  about  the  discrete  solution,  uh-  Then,  to  minimize  the  error 
due  to  this  approximation,  is  set  to  the  discrete  adjoint  solution  [40,  42].  That  is, 
^l^H  ^  is  chosen  such  that 

R%{\ih]^h,'^h)  =  R'h[^h]{^h,'^h)  -  J'[^h]{^h)  =  0,  VvH  G  Vh- 

Furthermore,  the  differences  u—uh  and  —  are  replaced  by  approximations  Uh  —  un 
and  —  'ipH^  where  and  i/j/j  are  surrogates  for  u  and  tp,  that  exist  in  a  richer  space, 
Vh,  than  Uh  and  ipn-  ti^is  work,  if  Vh  is  given  by  the  space  of  polynomials  of  order  p  in 
reference  space, 

Vh  =  {v  G  [L\n)Y  I  V  O  /,  G  [P^Kref)Y,  Vk  G  Th}  , 
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then  Vh  is  taken  to  be  the  space  of  polynomials  of  order  p  +  1.  In  particular, 


V/,  =  {V  G  [L\n)Y  I  VO  /,  G  [PP+\Kref)]\  Vk  G  Th]  . 

For  linear  geometry  elements  (i.e.  is  affine),  one  common  approach  is  to  define  u/j 
and  ■j/’/i  by  ^  local  or  patch  reconstruction  [72,  42].  For  example,  the  patch 
reconstruction  takes  the  following  form:  for  each  element  k,  find  v  G  such  that 

V  =  argmin  jjv  - 

where  V{k)  denotes  the  union  of  n  and  the  elements  that  share  a  face  with  k,  as  depicted 
in  Figure  4-1.  Then,  the  reconstruction  on  k  is  defined  by  restricting  v  to  k:  u/jj^  =  vj^. 


Figure  4-1:  The  patch  of  elements  centered  at  k,  denoted  'P(k),  highlighted  in  blue 

However,  in  the  curved  element  case,  since  pth  order  polynomials  in  the  reference  space 
are  generally  not  pth  order  polynomials  in  physical  space,  it  is  not  likely  that  u/j  dehned 
in  this  manner  is  in  the  space  Vh-  Furthermore,  the  patch  reconstruction  fails  to  account 
for  the  physics  of  the  problem.  Specifically,  the  reconstruction  does  not  properly  account 
for  the  propagation  of  information  in  convection  dominated  problems.  Thus,  while  it  is 
possible  to  generalize  the  patch  reconstruction  procedure  to  generate  a  reconstruction  in 
the  desired  space  or  to  chose  a  different  space  Vh,  a  different  technique  is  applied  here. 

To  motivate  the  procedure,  one  would  expect  that  choosing  u/j  and  xph  to  be  the  DG 
solutions  of  the  primal  and  adjoint  problems  posed  on  the  solution  space  Vh  would  give 
accurate  error  estimates.  However,  for  the  purposes  of  error  estimation,  these  solutions 
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are  regarded  as  too  expensive  to  compute.  Thus,  instead  of  computing  these  solutions,  the 
exact  solution  surrogates  are  defined  by  injecting  uh  and  into  Vh  and  taking  a  small 
number  of  element-block  Jacobi  iterations  on  the  p  +  1  discrete  problem.  To  be  precise,  the 
algorithm  has  the  following  steps: 

1.  Compute  the  injections  of  uh  and  into  Vh-  That  is,  find  G  and  G 
such  that 

N  N 

u//(x)  =  Uf(t)i{yL),  =  Y  W’ 

i=l  i=l 

where  {</>»}  for  i  =  . . .  ,N  forms  a  basis  of  Vh- 

2.  Take  K  element-block  Jacobi  iterations  on  the  primal  problem.  That  is,  for  k  = 

0,...,K-1, 

where  R{U^)  is  the  primal  residual  vector  evaluated  at  and  P([/^)  is  the  element- 
block  diagonal  of  the  Jacobian  matrix,  evaluated  at  - 

3.  Set  u,i(x)  =  Ya=i 

4.  Take  K  element-block  Jacobi  iterations  on  the  dual  problem.  That  is,  for  k  = 

0,...,K-1, 

q,k+i  ^q,k  _ 

where  R^{U^ is  the  residual  vector  for  the  dual  problem  taken  about  and 
evaluated  at 

5.  Set  'iphi^)  =  Eili 

Finally,  given  the  surrogate  solutions  Uh  and  ijjh,  the  error  can  be  approximated  using 
either  the  primal  residual  or  the  dual  residual; 

\J{u)-J{uh)\  ^  ^prim  —  \Rh{uH,'4’h-'4’H)\,  (4-5) 

iJ'(u)  -  J'(u/^)|  edual^\Rt{^H;'Vih-UH,'lpH)\7  (4.6) 

where  Rh  and  i?)(’  denote  the  primal  and  dual  residual  forms,  respectively,  on  the  space  Vh- 

An  elemental  error  indicator  is  required  for  the  mesh  adaptation  procedure.  This  indi¬ 
cator  is  constructed  by  averaging  the  primal  residual  and  dual  residual  approximations  to 
the  error  contribution  from  the  element.  In  particular,  for  each  element  k, 

^  (|i?h(u//,  i'lph  -  ^//)U)|  +  {uh  -  .  (4.7) 
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Another  error  estimate  can  be  constructed  by  summing  the  elemental  contributions  above: 


e  = 

K. 

Due  to  the  absolute  values,  this  estimate  is  generally  larger  than  both  €prim  and  e^uai- 

4.3  Numerical  Results 

In  this  section,  the  error  estimation  algorithm  is  applied  to  three  high  Re  test  cases.  In  all 
cases,  uniform  refinement  studies  have  been  conducted  to  enable  comparison  of  the  error 
estimates  to  the  error  computed  relative  to  fine  mesh  or  extrapolated  results.  Only  results 
for  the  standard  weighting  and  dual  consistent  discretizations  are  presented  because  the 
mixed  formulation  implementation  has  not  been  parallelized. 

4.3.1  Flat  Plate 

The  hrst  test  case  is  Rcc  =  1  x  lO'^,  M^o  =  0.25  flow  over  a  flat  plate  with  zero  pressure 
gradient.  The  coarse  mesh  in  this  study  is  a  1516  element  mesh  generated  by  mesh  adap¬ 
tation  using  the  p  =  3,  dual  consistent  discretization.  At  each  refinement  level,  the  mesh 
spacing  is  halved,  leading  to  a  set  of  three  nested  meshes.  These  meshes — referred  to  as  the 
coarse,  medium,  and  fine  meshes — contain  1516,  6064,  and  24256  elements.  The  exact  drag 
is  taken  to  be  that  computed  from  the  solution  of  the  p  =  4,  dual  consistent  discretization 
on  a  uniform  refinement  of  the  fine  mesh.  This  value  is  28.579496056  counts. 

Dual  Consistent  Discretization 

Figure  4-2  shows  the  drag  error  for  the  dual  consistent  discretization.  The  discretization 
attains  optimal  order  of  accuracy  for  the  drag  for  p  =  1,2.  However,  for  p  =  3,  the  order  is 
suboptimal.  The  cause  of  this  behavior  is  uncertain,  but  it  may  be  linked  to  the  singularity 
in  the  solution  at  the  leading  edge  of  the  plate  or  to  a  lack  of  resolution  at  the  boundary 
layer  edge.  Despite  the  observed  rate,  the  accuracy  of  the  p  =  3  solution  is  very  good. 
In  particular,  the  error  on  the  coarse  mesh  is  approximately  1  x  10“^%  of  the  total  drag. 
Figure  4-3  shows  the  error  estimates  Cprim  and  e^uai-  The  trends  observed  in  the  error 
estimates  are  very  similar  to  those  of  the  exact  error.  Moreover,  taking  the  maximum  of 
the  primal  and  dual  estimates  gives  a  slightly  conservative  but  still  accurate  error  estimate. 
Figure  4-4  shows  this  maximum  estimate,  as  well  as  its  effectivity,  defined  as  the  ratio  of 
the  maximum  estimate  to  the  actual  error.  The  effectivity  in  this  case  is  in  the  range  (1.5, 
4.6). 
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Figure  4-2:  Drag  error  for  the  dual  consistent  discretization  versus  uniform  grid  refinement 
for  the  flat  plate  test  case 


(a)  Primal  Residual  Estimate  (b)  Dual  Residual  Estimate 


Figure  4-3:  Primal  and  dual  residual  drag  error  estimates  for  the  dual  consistent  discretiza¬ 
tion  versus  uniform  grid  refinement  for  the  flat  plate  test  case 
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(a)  Maximum  Estimate 


(b)  Effect ivity 


Figure  4-4:  Maximum  error  estimate  and  effectivity  for  the  dual  consistent  discretization 
versus  uniform  grid  refinement  for  the  flat  plate  test  case 


Standard  Weighting  Discretization 

Figure  4-5  shows  the  drag  error  for  the  standard  weighting  discretization.  The  discretization 
attains  optimal  order  of  accuracy  for  the  drag  for  p  =  1.  For  p  =  2,  the  discretization  is 
suboptimal.  This  behavior  is  expected  and  agrees  with  the  results  observed  for  the  model 
problem  in  Chapter  3.  For  p  =  3,  the  order  is  also  suboptimal.  As  with  the  p  =  3,  dual 
consistent  discretization,  the  cause  of  this  behavior  is  unclear.  Figure  4-6  shows  the  error 
estimates  eprim  and  eduai-  Unlike  the  dual  consistent  discretization,  the  dual  inconsistent 
scheme  maximum  error  estimates  do  not  provide  a  conservative  estimate  for  all  p.  In 
particular,  as  shown  in  Figure  4-7,  the  effectivity  for  p  =  2  is  less  than  one  for  each  grid, 
and  it  decreases  with  grid  refinement. 


4.3.2  NACA  0012 

The  second  test  case  is  Rcc  =  1  x  10^,  Moo  =  0.25,  a  =  0  flow  over  a  NACA  0012  airfoil. 
The  computational  domain  and  meshes  for  this  case  are  as  described  in  Section  3.4.2.  The 
“exact”  drag  error  results  for  this  case  are  also  shown  in  Section  3.4.2.  Thus,  only  the  error 
estimates  and  effectivities  are  shown  here. 
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Figure  4-5:  Drag  error  for  the  standard  weighting  discretization  versus  uniform  grid  refine¬ 
ment  for  the  flat  plate  test  case 


(a)  Primal  Residual  Estimate  (b)  Dual  Residual  Estimate 


Figure  4-6:  Primal  and  dual  residual  drag  error  estimates  for  the  standard  weighting  dis¬ 
cretization  versus  uniform  grid  refinement  for  the  flat  plate  test  case 
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max(e  .  ,  e_,  ,)  (counts) 

'  pnm  dual'  '  ' 


(a)  Maximum  Estimate  (b)  Effectivity 


(c)  Effectivity  (zoom) 


Figure  4-7:  Maximum  error  estimate  and  effectivity  for  the  standard  weighting  discretization 
versus  uniform  grid  refinement  for  the  flat  plate  test  case 
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Dual  Consistent  Discretization 


Figure  4-8  shows  the  primal  and  dual  error  estimates.  Unlike  the  flat  plate  results,  the  two 
estimates  show  somewhat  different  behavior.  For  example,  on  the  second  refinement,  the 
p  =  2  primal  estimate  decreases  dramatically  while  the  dual  estimate  is  relatively  constant. 
However,  as  Figure  4-9  shows,  the  maximum  of  the  two  error  estimates  generally  provides  a 


^  °  logjiho/h) 


(a)  Primal  Residual  Estimate 


(b)  Dual  Residual  Estimate 


Figure  4-8;  Primal  and  dual  residual  drag  error  estimates  for  the  dual  consistent  discretiza¬ 
tion  versus  uniform  grid  refinement  for  the  NACA  0012  test  case 

conservative  but  reasonably  accurate  estimate.  In  particular,  the  effectivity  of  the  maximum 
estimate  is  between  0.8  and  2.1  for  all  cases  except  p  =  2  on  the  medium  mesh. 

In  Section  3.4.2,  it  was  observed  that  the  order  of  accuracy  for  p  >  3  for  this  case  is  lower 
than  that  expected  for  smooth  solutions.  This  behavior  was  attributed  to  under-resolution 
at  the  boundary  layer  edge.  To  examine  this  claim  in  more  detail.  Figures  4-10  and  4-11 
show  the  elemental  error  estimate,  and  Mach  number  for  p  =  2  and  p  =  4,  respectively. 
The  figure  demonstrates  that,  in  the  p  =  2  case,  while  the  error  estimate  is  non-zero  at  the 
boundary  layer  edge,  the  errors  are  largest  in  the  flow  upstream  of  the  airfoil.  Alternatively, 
in  the  p  =  4  case,  the  error  estimate  is  largest  near  the  boundary  layer  edge. 

Standard  Weighting  Discretization 

Figure  4-12  shows  the  primal  and  dual  error  estimates.  The  results  are  qualitatively  similar 
to  those  for  the  dual  consistent  discretization.  The  maximum  error  estimate  and  effectivity 
are  shown  in  Figure  4-13.  As  with  the  dual  consistent  discretization,  the  maximum  error 
estimate  generally  provides  an  accurate  error  estimate  for  this  case.  However,  for  p  =  2 
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max(e  .  ,,  e_,  ,)  (counts) 
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(a)  Maximum  Estimate 


(b)  Effectivity 


(c)  Effectivity  (zoom) 


Figure  4-9:  Maximum  error  estimate  and  effectivity  for  the  dual  consistent  discretization 
versus  uniform  grid  refinement  for  the  NACA  0012  test  case 
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Figure  4-10:  Elemental  error  estimate  (left)  and  Mach  number  (right)  for  the  p  =  2,  dual 
consistent  discretization  on  the  coarse  mesh  for  the  NACA  0012  test  case 
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Figure  4-11:  Elemental  error  estimate  (left)  and  Mach  number  (right)  for  the  p  =  4,  dual 
consistent  discretization  on  the  coarse  mesh  for  the  NACA  0012  test  case 
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(a)  Primal  Residual  Estimate 


(b)  Dual  Residual  Estimate 


Figure  4-12:  Primal  and  dual  residual  drag  error  estimates  for  the  standard  weighting 
discretization  versus  uniform  grid  refinement  for  the  NACA  0012  test  case 

the  error  estimate  is  converging  more  rapidly  than  the  actual  error,  leading  to  a  decreasing 
effectivity.  Thus,  on  the  fine  mesh,  the  p  =  2  effectivity  is  approximately  0.5.  The  only 
other  case  with  an  effectivity  outside  the  range  of  (0.8,  2)  is  the  p  =  3  fine  grid  result,  where 
the  effectivity  is  almost  seven. 

4.3.3  RAE  2822 

The  final  test  case  is  Rcc  =  6.3  x  10®,  Mqo  =  0.6,  a  =  2.57°  flow  over  a  RAE  2822  airfoil.  The 
computational  domain,  meshes,  and  “exact”  drag  for  this  case  are  given  in  Section  3.4.3, 
which  also  contains  the  “exact”  drag  error  results.  Thus,  only  the  error  estimates  and 
effectivities  are  shown  here.  Furthermore,  only  the  dual  consistent  discretization  is  used. 

Figure  4-14  shows  the  primal  and  dual  error  estimates.  The  primal  error  estimate  for 
p  =  1  on  the  coarse  mesh  is  so  large  (3.5  x  10^  counts)  that  it  is  not  shown  in  the  figure.  This 
large  error  estimate  results  from  the  fact  that  the  element-block  Jacobi  iterative  procedure 
used  to  compute  the  surrogate  primal  solution  diverges  for  this  case.  However,  for  the 
other  cases,  where  the  element-block  Jacobi  solve  is  well-behaved,  the  error  estimates  are 
reasonable.  The  maximum  error  estimate  and  its  effectivity  are  shown  in  Figure  4-15.  The 
effectivity  is  in  the  range  (0.9, 11.0)  for  all  but  the  p  =  1  result  on  the  coarse  mesh. 

As  noted  in  Section  3.4.3,  the  order  of  accuracy  for  p  >  3  for  this  case  is  below  that 
expected  for  smooth  flow.  Figure  4-16  shows  the  p  =  3  elemental  error  estimates  on  the 
coarse  and  fine  meshes.  At  first  glance,  it  appears  that  the  error  is  small  everywhere  near 
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(a)  Maximum  Estimate 


(b)  Effectivity 


(c)  Effectivity  (zoom) 


Figure  4-13:  Maximum  error  estimate  and  effectivity  for  the  standard  weighting  discretiza¬ 
tion  versus  uniform  grid  refinement  for  the  NACA  0012  test  case 
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(a)  Primal  Residual  Estimate 


(b)  Dual  Residual  Estimate 


Figure  4-14:  Primal  and  dual  residual  drag  error  estimates  for  the  dual  consistent  discretiza¬ 
tion  versus  uniform  grid  refinement  for  the  RAE  2822  test  case 


(a)  Maximum  Estimate 


(b)  Effectivity 


Figure  4-15:  Maximum  error  estimate  and  effectivity  for  the  dual  consistent  discretization 
versus  uniform  grid  refinement  for  the  RAE  2822  test  case 
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the  airfoil.  However,  close  inspection  of  the  trailing  edge  region,  as  shown  for  the  fine 
mesh  in  Figure  4-17,  shows  that  the  elemental  errors  in  this  region  are  large.  These  errors 
indicate  a  lack  of  resolution  in  this  region.  Thus,  one  can  conclude  that  under-resolution  of 
the  trailing  edge  singularity  is  leading  to  a  loss  of  accuracy  for  p  >  3  for  this  case. 
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(a)  Coarse  mesh 


(b)  Fine  mesh 


Figure  4-16:  Elemental  error  estimate  for  the  p  =  3,  dual  consistent  discretization  on  the 
coarse  and  fine  meshes  for  the  RAE  2822  test  case 
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(a)  Trailing  edge 


(b)  Trailing  edge  zoom 


Eigure  4-17:  Elemental  error  estimate  at  the  trailing  edge  for  the  p  =  3,  dual  consistent 
discretization  on  the  fine  mesh  for  the  RAE  2822  test  case 
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Chapter  5 


Output-Based  Adaptation 


In  this  chapter,  the  output  error  estimate  is  used  to  drive  a  mesh  adaptation  algorithm. 
The  algorithm  combines  the  output  error  estimate  with  interpolation  error  estimates  to  de¬ 
termine  a  desired  metric  for  the  new  mesh.  This  technique,  described  in  Section  5.1,  closely 
follows  those  of  Venditti  and  Darmofal  [102,  103]  and  Fidkowski  and  Darmofal  [42,  39]. 
Thus,  the  main  contributions  of  this  chapter  are  the  application  of  the  algorithm  to  high- 
order  discretizations  of  the  RANS  equations  and  the  use  of  curved,  boundary  conforming 
meshes.  Two  techniques  for  generating  curved  meshes  are  detailed  in  Section  5.1.3.  Sec¬ 
tion  5.2  provides  numerical  results. 

5.1  Algorithm 

The  adaptation  algorithm  used  here  takes  the  following  form: 

1.  Generate  an  initial  mesh. 

2.  Solve  the  steady  state  problem,  R(U)  =  0,  where  R  denotes  the  residual  vector  and 
U  denotes  the  discrete  solution.  Specifically, 

(a)  Set  the  solution  to  some  initial  guess  U'^.  Initialize  the  CFL  number,  and  com¬ 
pute  the  corresponding  time  step.  At.  Set  i  =  1. 

(b)  While  ||R(U*-i)||  Xltoi, 

i.  Update  the  primal  state:  U*  =  -f 

ii.  Update  At  according  to  Section  2.4. 
hi.  i  — >  i  -|-  1. 

where  6toi  is  the  solution  tolerance,  A4  is  the  mass  matrix. 


TU- 


R(U 
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0-p>  I  r\  j  I  'jp 

3.  Solve  the  adjoint  problem:  where  J  is  the  discrete  output  of  interest. 

4.  Using  U  and  'h,  compute  e,  the  output  error  estimate,  as  described  in  Chapter  4. 

5.  If  e  <  etoi,  where  etoi  is  the  requested  error  tolerance,  quit.  Otherwise,  adapt  the 
mesh,  and  return  to  step  2. 

Thus,  at  each  adaptation  iteration,  an  error  estimate  is  computed,  and,  assuming  that 
estimate  is  greater  than  the  desired  tolerance,  a  new  computational  mesh  is  generated.  The 
desired  mesh  is  specified  by  means  of  an  n  x  n  Riemannian  metric,  M,  defined  over  the 
physical  domain,  where  n  is  the  dimension  of  the  domain.  The  metric  contains  information 
on  the  desired  mesh  spacing  in  physical  space.  In  particular,  under  the  metric,  the  mesh 
spacing  is  required  to  be  unity  in  all  directions.  To  be  precise,  let  t  be  a  unit  vector,  and  let 
/i((x)  be  the  desired  mesh  spacing  in  the  t  direction  at  x.  Then,  the  metric  at  x  is  defined 
by  requiring 

=  1,  (5.1) 

for  all  directions  t.  This  requirement  defines  an  ellipsoid  in  physical  space.  The  eigenvectors 
of  M,  Cj  for  i  =  1, . .  .n,  give  the  directions  of  the  principal  axes  of  an  ellipsoid,  and  the 
eigenvalues  of  M,  A*  for  i  =  1, . . .  ,n  are  related  to  the  lengths  of  the  ellipsoid  axes  by 
h'f  =  l/\i. 

The  mesh  metric  calculation  involves  two  parts:  anisotropy  detection  and  mesh  opti¬ 
mization. 


5.1.1  Anisotropy  Detection 

The  desired  anisotropy  calculation  is  motivated  by  a  Hessian-based  analysis  [21,  103].  The 
analysis  begins  with  an  interpolation  error  estimate.  If  uh  is  a  linear  interpolant  of  a  smooth 
scalar  function  u,  assuming  the  interpolation  error  is  zero  at  the  segment  endpoints,  the 
interpolation  error  along  a  segment  dx  is  bounded  by 


max|n(x)  —  tt//(x)|  < 

xG5x  x€(5x 


U 


(2) 

5x 


X 


(2) 

where  /i^x  is  the  length  of  the  segment,  Ci  is  a  constant  that  is  independent  of  u,  and 
is  the  second  derivative  of  u  in  the  dx  direction  [96].  Using  the  Hessian  matrix,  H,  where 


d'^u 

dxidxj  ’ 
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the  interpolation  error  bound  can  be  rewritten  as 


max  |u(x)  —  ttH(x)|  <  Ci  max  ((5x'^|H(x)|(5x)  , 

xG5x  xGSx 

where  |H|  =  V|A|V-i  for  H  =  VAV-h 

The  desired  mesh  is  designed  to  give  equal  interpolation  errors  in  all  directions.  Thus, 
for  some  constant  C,  the  requested  mesh  spacing  should  satisfy 

h?F|H|t  =  i  (5.2) 

for  all  unit  vectors  t.  Using  (5.1),  (5.2)  can  be  satisfied  by  requiring  that 

M  =  C|H|.  (5.3) 


This  relationship  determines  the  metric  up  to  a  multiplicative  constant  which  is  determined 
by  the  mesh  optimization  step.  Furthermore,  given  that  the  exact  Hessian  matrix  is  not 
known  it  must  be  approximated  to  compute  the  desired  anisotropy. 

Following  Fidkowski  [42,  39],  this  analysis  can  be  extended  to  p  >  1  interpolation.  As 
with  the  p  =  1  case,  the  analysis  begins  with  an  interpolation  error  bound: 


max|u(x)  —  Uff(x)l  <  max 

xGSx  xGSx 


U 


(p+l) 

6x 


(5.4) 


where  uh  is  a  pth  order  interpolant  of  u  that  interpolates  u  exactly  at  the  endpoints  of 
the  segment  5x.  For  p  >  1,  there  is  no  analog  of  the  Hessian  matrix  that  can  be  used  to 
bound  the  interpolation  error  and  define  the  mesh  metric.  Instead,  the  principal  stretching 
directions  are  explicitly  computed.  Specifically,  define  ei  to  be  the  direction  of  the  maximum 
(p  +  1)  derivative.  Let  e2  be  the  direction  of  the  maximum  (p  +  1)  derivative  in  the  plane 
orthogonal  to  ei,  and  so  on.  Then,  the  desired  anisotropy  is  determined  by  requiring  equal 
interpolation  error  in  the  directions  given  by  e*.  Thus,  for  i,  j  =  1, . . . ,  n. 


C  ^ 


hj 


u^\ 

4r‘V 


i/(p+i) 


As  with  (5.3),  this  relationship  determines  the  metric  up  to  a  multiplicative  constant.  This 
constant  is  based  on  the  error  estimate,  as  shown  in  Section  5.1.2.  Furthermore,  as  with  the 
exact  Hessian,  the  exact  (p  +  1)  derivatives  are  not  known.  Thus,  they  must  be  estimated 
to  compute  the  desired  anisotropy. 

To  apply  this  anisotropy  detection  method  on  curved  meshes,  a  minor  modification 
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is  required.  Specifically,  given  that  the  numerical  solution  is  a  pth  order  polynomial  in 
reference  space  but  not  physical  space,  the  interpolation  error  estimate  (5.4)  does  not  hold. 
However,  an  analogous  bound  holds  in  reference  space.  In  particular,  let 

u(x)  =  (uo/)(^)  =  u(^), 


where  /  is  the  map  from  reference  to  physical  space  and  ^  denotes  the  position  in  the 
reference  space.  Then,  if  uh  is  a  pth  order  interpolant  of  u. 


max  |u(^)  -  <  C2h^^t^  max 


u 


(p+i) 


where  6^  is  a  segment  in  reference  space.  Thus,  following  the  technique  above,  one  can  define 
the  principal  stretching  directions  in  reference  space,  e*.  Then,  requiring  equal  interpolation 
errors  in  the  principal  directions  gives 

h  ■  ( 

=  constant  ^  =  \  , 

hj  14^  ^ 


j  i/(p+i) 


where  hi  are  the  desired  mesh  spacings  in  reference  space.  This  relationship  determines  the 
metric  in  reference  space  to  within  a  multiplicative  constant. 


To  define  the  metric  in  physical  space,  the  mapping  from  the  reference  to  the  physical 
space  is  used.  In  particular,  given  the  metric  in  reference  space,  M,  at  ^q,  the  points  of  the 
corresponding  ellipsoid,  ^g,  satisfy 

(C-^o)'^M(^0)(^e-^o)  =  l. 

The  image  of  this  ellipsoid  under  the  map  /  is  given  by  /(^g).  Certainly,  since  /  may  be 
nonlinear,  this  image  may  not  be  an  ellipsoid.  However,  linearizing  /  about  gives 

/(O^Xg  =  /(Co)+I)/(^o)(^e-^o)- 

Using  this  linearized  form,  the  points  Xg  define  an  ellipsoid  centered  at  xq  =  /(^o)-  This 
ellipsoid  is  used  to  define  the  metric  in  physical  space  at  the  point  xq.  To  illustrate  this 
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calculation,  let 


hi 

H  =  VAV"^  = 

ei 

hn 

- 1 

iq; 

_ 1 

If  t  is  a  unit  vector,  then 

^e-^0  =  Ht, 

gives  the  point  on  the  ellipsoid  defined  by  the  metric  in  reference  space.  Then, 

xg  -  /(^o)  =  ^/(^o)Ht, 

gives  the  ellipsoid  defining  the  metric  in  physical  space.  Thus,  the  directions  and  lengths  of 
the  principal  axes  of  this  ellipsoid  are  given  by  the  left  singular  vectors  and  singular  values 
of  ZI/(^q)H.  This  process  is  illustrated  in  two  dimensions  in  Figure  5-1. 


Reference  Space  Physical  Space 


Figure  5-1:  Illustration  of  transfer  of  desired  mesh  anisotropy  from  reference  to  physical 
space  for  two  dimensional  case 

As  noted  earlier,  the  exact  {p  +  1)  derivatives  are  not  known.  Moreover,  given  that  the 
discrete  solution  uh  is  order  p,  it  contains  no  information  on  the  {p+  1)  derivatives.  Thus, 
in  this  work,  the  {p  +  1)  derivatives  are  estimated  using  the  surrogate  solution,  u/^,  which 
is  computed  for  the  error  estimate.  This  function  is  a  (p  -|-  l)th  order  polynomial  in  the 
reference  element.  Thus,  the  {p  +  1)  derivatives  are  constant  over  the  reference  element, 
which  implies  that  the  desired  anisotropy  in  the  reference  element  is  also  constant.  To 
transfer  the  anisotropy  request  to  physical  space,  the  Jacobian,  Df,  is  required.  In  general, 
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this  Jacobian  varies  over  the  reference  element.  In  this  work,  the  value  at  the  centroid  of 
the  reference  element  is  used  to  calculate  the  desired  anisotropy  in  the  physical  space. 


5.1.2  Mesh  Optimization 


As  noted  in  Section  5.1.1,  the  anisotropy  detection  step  only  determines  the  desired  mesh 
metric  to  within  a  multiplicative  constant.  That  constant  is  determined  by  the  mesh  op¬ 
timization  step.  Mesh  optimization  refers  to  determining  which  elements  to  refine/coarsen 
and  the  appropriate  level  of  refinement/coarsening.  In  this  work,  the  method  of  Fid- 
kowski  [42,  39]  is  used. 

The  main  idea  of  the  method  is  to  attempt  to  equidistribute  the  error  over  the  desired 
mesh.  To  make  this  idea  concrete,  let  eo  be  the  desired  global  error.  This  desired  error  is 
given  by 

eo  =  max{r]a€,r]t€toi), 

where  e  =  is  the  conservative  error  estimate  for  the  solution  on  the  current  mesh, 

etoi  is  the  error  tolerance,  and  0  <  ?7a  <  1  and  0  <  <  1  are  the  adaptation  aggressiveness 

and  target  error  fraction,  respectively.  Then,  letting  be  the  estimated  total  number  of 
elements  in  the  desired  mesh,  eo/Nf  is  the  allowable  error  per  element  in  the  desired  mesh. 

Furthermore,  let  be  the  number  of  elements  in  the  desired  mesh  contained  in  the 
element  k  in  the  current  mesh.  Then,  can  be  approximated  by 


riK 


2=1 


(5.5) 


where  hi  denotes  the  desired  mesh  spacings  and  h'^  denotes  the  current  mesh  spacings.  The 
current  mesh  spacings  are  computed  from  the  current  grid-implied  metric.  To  define  the 
grid-implied  metric,  let  i  be  the  affine  map  from  the  reference  element  to  the  unit  equi¬ 
lateral  triangle/tetrahedron.  This  map,  together  with  the  map  from  the  reference  element 
to  physical  space,  is  illustrated  in  Figure  5-2.  The  grid-implied  metric  is  dehned  by  the 
left  singular  vectors  and  singular  values  of  the  product  where  Df^  is  evaluated 

at  the  centroid  of  the  reference  element.  This  definition  is  an  extension  of  that  used  by 
Fidkowski  [42,  39]  to  the  curved  element  case. 

Given  n^,  one  can  write  the  allowable  error  for  the  elements  in  the  desired  mesh  con¬ 
tained  in  the  element  n  of  the  current  mesh  as  nf^eo/Nf.  To  set  the  mesh  spacing  on  the 
desired  mesh  to  achieve  this  error,  it  is  necessary  to  relate  the  mesh  spacing  to  the  expected 
error.  This  relationship  is  given  by  an  a  priori  error  estimate.  Specifically,  it  is  assumed 
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Unit  Equilateral  „  ^  t-  •  i  nu  •  i  ci 

^  Reierence  Triangle  Physical  Element 


Figure  5-2:  Mappings  from  the  reference  element  to  the  unit  equilateral  triangle  and  physical 
space 


that 

WJ  ’ 

where  e%  is  the  current  error  estimate  on  k,  given  by  4.7,  and  is  the  expected  error  on 
the  elements  of  the  desired  mesh  contained  in  k.  In  this  estimate,  it  is  assumed  that  the 
smallest  mesh  spacing  principal  directions  in  the  current  and  desired  meshes,  and  ei, 
align.  Clearly,  this  alignment  need  not  occur.  However,  one  expects  that,  as  the  adaptation 
progresses,  the  principal  directions  will  converge. 

The  order  of  accuracy  is  taken  to  be  =  s  +  t,  where  s  is  the  convergence  rate  of  the 
primal  solution  in  an  appropriate  energy  norm  and  t  is  the  convergence  rate  of  the  adjoint 
solution  in  an  appropriate  energy  norm.  See  Appendix  C  for  an  a  priori  output  error  esti¬ 
mate  motivating  this  choice.  For  elliptic  problems,  assuming  exact  geometry  representation 
and  a  solution  space  containing  pth  order  polynomials  in  physical  space,  one  can  show 

s  =  min(p,  7  —  1),  t  =  min(p,  —  1), 


where  7  and  7^  denote  the  regularities  of  the  exact  primal  and  dual  solutions,  respectively — 
i.e.  u  G  and  'ip  G  (U).  In  this  work,  unless  otherwise  noted,  optimal  rates — 

s  =  t  =  p — are  assumed  for  all  elements  not  containing  geometric  singularities.  For  elements 
touching  a  geometric  singularity — e.g.  the  leading  edge  of  an  infinitely  thin  flat  plate  or  a 
sharp  trailing  edge — s  -|-  t  =  1  is  used. 


Equating  the  allowable  and  expected  errors  on  the  element  k  gives 


ep 

'Nf 


= 


(5.6) 


Using  (5.5),  hi /hi  can  be  written  in  terms  of  the  current  mesh  spacing  ratio,  the  desired 
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mesh  spacing  ratio  from  the  anisotropy  detection  step,  and  n^.  In  two  dimensions, 


hi  /  1  /i2  hi  \ 

hi  \n^hlh2) 

Then,  substituting  (5.7)  into  (5.6)  and  solving  for  gives 

n  -  rA 

"  Uo  V  [hlh2) 

Finally,  summing  over  the  elements  gives  an  equation  for  Nf. 

kGTh  K.&Th 


eo 


hi  hi 
hi  h2 


T'ft+2 


(5.7) 


(5.8) 


(5.9) 


Equation  (5.9)  can  be  solved  for  Nf,  allowing  to  be  computed  from  (5.8).  Then,  hi 
can  be  computed  using  (5.7).  This  calculation  combined  with  the  desired  anisotropy  from 
Section  5.1.1  completely  determines  the  desired  mesh  metric. 


5.1.3  Curved  Mesh  Generation 

The  metric  computed  as  shown  in  Sections  5.1.1  and  5.1.2  can  be  input  to  a  mesh  generation 
package  to  generate  a  linear  mesh.  In  this  work,  the  Bi-dimensional  Anisotropic  Mesh 
Generator  (BAMG)  package  is  used  [19,  58].  However,  it  is  well  known  that  to  fully  realize 
the  potential  of  high-order  methods,  the  geometry  representation  should  be  high-order 
accurate  as  well  [13].  Thus,  for  curved  geometries,  linear  meshes  are  not  appropriate. 
Unfortunately,  curved  mesh  generation  is  a  non-trivial  task.  In  this  work,  curved  meshes 
are  generated  by  first  generating  a  linear  mesh  and  then  operating  on  that  linear  mesh  to 
arrive  at  a  suitable  curved  mesh.  Specifically,  two  methods  will  be  examined:  user-generated 
global  mappings  and  linear  elasticity  mesh  movement. 


User-Generated  Global  Mapping 

The  user-generated  global  mapping  approach  is  motivated  by  the  observation  that  it  is 
often  possible  to  specify  a  global  map  from  a  simple,  linear  meshing  domain — such  as  a 
rectangle — to  the  physical  domain  of  interest.  Figure  5-3  shows  a  cartoon  example  of  such 
a  map,  g  :  VLm  where  denotes  the  meshing  domain. 

Given  this  map,  the  idea  applied  here  is  to  generate  a  linear  mesh  of  the  meshing  domain 
and  then  map  that  mesh  to  the  physical  domain.  Of  course,  the  resulting  mesh  is  required 
to  satisfy  the  desired  mesh  metric  in  physical  space.  To  enforce  this  condition,  a  metric  is 
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Figure  5-3:  Illustration  of  the  global  map,  g,  from  meshing  to  physical  domain 

generated  for  the  meshing  domain.  The  procedure  for  transferring  the  metric  in  physical 
space  to  a  metric  in  meshing  space  is  analogous  to  that  of  transferring  the  anisotropy  request 
from  the  reference  to  the  physical  space  as  shown  in  Section  5.1.1.  Specifically,  ignoring  any 
high-order  terms,  the  map  g~^  takes  the  ellipsoid  defined  by  the  metric  in  physical  space 
to  an  ellipsoid  in  meshing  space: 

Fe  =  +  T>5f"^(xo)(Xe  “  Xq), 

where  xq  is  the  center  of  the  ellipsoid  in  physical  space,  Xg  are  points  on  the  ellipsoid  in 
physical  space,  and  re  are  points  on  the  ellipsoid  in  meshing  space.  This  ellipsoid  is  used 
to  define  the  mesh  metric  in  meshing  space. 

Using  this  metric,  one  can  generate  a  linear  mesh  of  the  meshing  space.  Then,  a  high- 
order  mesh  of  physical  space  is  generated  by  adding  uniformly  spaced  higher-order  nodes 
to  the  linear  mesh  in  meshing  space  and  mapping  the  nodes  to  the  physical  domain.  After 
computing  the  nodes  in  the  physical  domain,  the  mesh  of  the  physical  domain  is  defined  by 
high-order  Lagrange  interpolation  of  these  nodes.  This  algorithm  restricts  the  use  of  the 
map  g  and  the  meshing  domain  to  the  mesh  generation  step.  The  alternative  is  to  define 
the  mesh  of  the  physical  space  as  the  image  of  the  linear  mesh  of  the  meshing  domain 
under  the  map  g.  While  this  latter  approach  has  some  nice  features — e.g.  exact  geometry 
representation — it  is  not  used  here  because  it  requires  the  map  g  to  compute  the  residual. 

Linear  Elasticity  Node  Movement 

While  the  user-generated  mapping  technique  yields  good  results,  it  places  a  fairly  high 
burden  on  the  user  to  generate  a  suitable  meshing  domain  and  associated  global  map. 
Especially  for  complicated  three-dimensional  problems,  this  burden  is  unacceptable.  Thus, 
a  second,  more  general,  approach  is  also  explored.  In  this  approach,  the  metric  is  used 
directly  to  generate  a  linear  mesh  of  the  physical  domain.  After  adding  uniformly  spaced 
higher-order  nodes  to  the  linear  mesh,  all  mesh  nodes  are  moved  to  generate  a  valid,  high- 
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order  mesh. 

The  node  movement  is  governed  by  the  linear  elasticity  equations,  as  used  by  [77].  Let 
Ql  denote  the  domain  covered  by  the  linear  triangulation,  T/j  /,.  Then,  in  two  dimensions, 
the  node  movement  is  governed  by 

V  •  (VV  +  <t)  =  0,  X  G  rii 
V-Vsc  =  0,  xgSLIl. 

where  V  =  [u,  v]'^  is  the  displacement  vector, 

AR{V  •  V)  0 
0  AR{V  •  V) 

and  AR  is  the  aspect  ratio  of  the  underlying  linear  element.  Using  the  aspect  ratio  to  set 
the  material  properties  in  this  manner  results  in  high  aspect  ratio  cells  that  act  like  nearly 
incompressible  materials,  which  adds  robustness  to  the  node  movement  in  highly  stretched 
regions  [77].  The  Dirichlet  boundary  conditions,  V^C;  are  set  by  requiring  that 

Xp  =  X  +  Vsc'(x),  X  G 

where  Xp  is  the  projection  of  the  point  x  onto  the  exact  geometry  of  the  physical  domain. 

The  node  movement  equations  are  discretized  using  a  high-order  continuous  Galerkin 
discretization.  The  order  of  this  discretization,  denoted  q,  is  set  by  the  desired  order  of  the 
geometry  representation.  Thus,  the  mesh  movement  is  given  by  the  solution  of  the  following 
problem:  find  V/j  G  X'^  such  that 


where  the  function  spaces  are  given  by 

Xl  =  {cl>h  G  C\VLl)  I  Mn  G  P^{n),  Vk  G  Th,L}. 

Xl^  =  {</>/,  G  I  =  0}, 

Kbc  -  {^h  G  Xl  I  =  Vbc(x,),  Vj  =  1, . . . ,  Xfe,}, 

and  {xj  I  j  =  l,...,Nbn}  denotes  the  nodes  on  dQi-  Clearly,  this  discrete  problem  is 
linear.  Thus,  computing  the  node  movement  solution  requires  a  single  linear  solve.  In  this 
implementation,  the  sparse  linear  system  is  solved  using  Gaussian  elimination. 
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5.2  Numerical  Results 


In  this  section,  the  adaptation  algorithm  is  applied  to  four  two-dimensional  test  flows:  a  flat 
plate,  an  ellipse,  and  two  airfoil  flows.  In  all  cases,  the  output  used  to  drive  the  adaptation 
is  drag,  and  the  Mach  number  is  used  to  determine  the  anisotropy.  Furthermore,  the 
adaptation  aggressiveness  factor,  r]a,  is  set  to  O.I,  and  the  error  tolerance  is  set  arbitrarily 
small  such  that  it  will  never  be  reached.  This  tolerance  was  chosen  to  simulate  very  stringent 
accuracy  requirements.  For  brevity,  the  results  for  all  three  source  term  discretizations 
are  shown  only  for  the  flat  plate  case.  For  the  remaining  cases,  only  the  dual  consistent 
discretization  is  used. 

5.2.1  Flat  Plate 

The  hrst  test  case  is  M^o  =  0.25,  Rcc  =  1  x  10^  flow  over  a  flat  plate  with  zero  pressure 
gradient.  Figure  5-4  shows  the  adaptation  histories  for  each  source  term  discretization 
investigated.  In  particular,  the  error  estimate  and  actual  error  are  plotted  versus  DOF  for 
each  adaptation  iteration.  The  error  estimate  shown  is  the  maximum  of  the  primal  and 
dual  residual  estimates  given  by  (4.5)  and  (4.6),  respectively.  The  actual  error  is  computed 
relative  to  the  same  value  used  for  this  case  in  Section  4.3.1. 

From  the  figure,  it  is  clear  that,  for  all  the  source  term  discretizations  examined,  the 
p  =  2  and  p  =  3  discretizations  are  substantially  more  efficient  than  p  =  1  in  terms  of  DOF 
required  to  achieve  a  desired  accuracy.  For  example,  using  the  dual  consistent  discretization, 
the  p  =  2  and  p  =  3  discretizations  both  achieve  an  error  of  approximately  5  x  10“^  drag 
counts  using  approximately  2.5  x  10^  DOF.  Alternatively,  the  p  =  I  discretization  requires 
more  than  4  x  10^  DOF,  an  order  of  magnitude  increase. 

To  facilitate  comparison  of  the  three  source  discretizations,  the  adaptation  histories  are 
shown  again  in  Figure  5-5.  For  all  p,  the  three  source  term  discretizations  produce  very 
similar  error  estimates.  However,  the  actual  errors  are  somewhat  different.  The  most  drastic 
differences  occur  for  p  =  2,  where  the  dual  consistent  and  mixed  formulation  discretizations 
significantly  outperform  the  standard  weighting  scheme.  For  example,  using  7374  DOF,  the 
standard  weighting  scheme  gives  an  error  of  approximately  4  x  10“^  counts  while  the  dual 
consistent  and  mixed  formulations  give  errors  of  less  than  3  x  10“^  counts  using  similar 
DOF. 

Figures  5-6  through  5-8  show  skin  friction  distributions  on  the  initial  mesh  and  after 
two  adaptation  iterations  for  p  =  1,2,3.  The  figures  show  that  the  computed  skin  friction 
matches  the  experimental  data  of  Wieghart  [32]  very  closely  after  two  adaptation  iterations 
for  all  p  and  all  discretizations  used.  Furthermore,  the  skin  friction  distributions  reinforce 
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DOF 


(a)  Standard  weighting,  Error  estimate 


DOF 


(c)  Dual  consistent,  Error  estimate 


DOF 

(b)  Standard  weighting,  Actual  error 


(d)  Dual  consistent.  Actual  error 


(e)  Mixed  formulation.  Error  estimate  (f)  Mixed  formulation,  Actual  error 


Figure  5-4:  Estimated  and  actual  drag  error  versus  DOF  for  adaptation  on  the  drag  on  a 
flat  plate,  comparison  of  p  =  1,  2,  3  results 
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DOF  DOF 


(a)  p  =  1,  Error  estimate 


(b)  p  =  1,  Actual  error 


(c)  p  =  2,  Error  estimate 


(d)  p  =  2,  Actual  error 


(e)  p  =  3,  Error  estimate 


(f)  p  =  3,  Actual  error 


Figure  5-5:  Estimated  and  actual  drag  error  versus  DOF  for  adaptation  on  the  drag  on  a  flat 
plate,  comparison  of  standard  weighting  (SW),  dual  consistent  (DC)  and  mixed  formulation 
(ME)  discretizations 
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the  conclusion  from  the  adaptation  histories  that  the  p  =  2  and  p  =  3  discretizations  are 
far  superior  to  p  =  1.  In  particular,  for  all  three  discretizations,  after  two  adaptation 
iterations,  the  p  =  2  and  p  =  3  skin  friction  distributions  appear  smoother  than  the  p  =  1 
result,  despite  the  fact  that  the  p  =  1  uses  more  than  twice  the  DOF. 

5.2.2  Ellipse 

The  second  test  case  is  Moo  =  0.25,  Rbc  =  1  x  10^,  0  =  0  flow  over  an  ellipse.  The  ellipse 
has  a  semi-major  axis  length  a  =  1.0  and  a  semi-minor  axis  length  b  =  0.1.  Furthermore, 
assuming  the  x-axis  is  oriented  along  the  major  axis  of  the  ellipse,  the  flow  is  symmetric 
about  y  =  0.  Thus,  only  the  upper  half  of  the  domain  is  used  for  the  simulation.  The 
boundary  of  the  physical  domain  is  given  by  dQ  =  U^^qFj,  where 


ho 

=  1 

-  100  <  X  <  -1,  y  =  0}, 

Fi 

=  {{x,y)  1 

—  l<x<l,  y  =  b\/l  —  x^/a^}. 

F2 

=  {ix,y)  1 

1  <  X  <  100,  y  =  0}, 

hs 

=  {ix,y)  1 

—  100  <  X  <  100,  y  =  V 100  —  x^} 

On  F3,  the  full  state  boundary  condition  is  used.  Boundaries  Fq  and  F2  are  symmetry  lines, 
and  Fi  is  a  no  slip,  adiabatic  wall. 

Given  the  simple  geometry  of  the  domain,  it  is  straightforward  to  construct  a  mapping 
from  a  rectangular  meshing  domain  to  the  physical  domain.  Specihcally,  the  mapping  from 
the  meshing  space  to  the  physical  space,  g  :  VLm  is  given  by 


exp,^  COST] 

9{C,V)  = 

exp^ 

logR 

byjl  (logii  C)-FCsm7?j 

where  R  =  100.0  and  Dm  =  I  0  <  ^  <  logii,  0  <  ??  <  tt}.  Furthermore,  superpara- 

metric  elements  are  used.  In  particular,  the  geometry  of  each  element  is  described  using 
q  =  5  Lagrange  polynomials  for  all  p.  To  illustrate  the  mesh  map.  Figure  5-9  shows  the 
initial  mesh  used  for  the  adaptation  run  in  both  the  meshing  and  physical  spaces. 

Figure  5-10  shows  the  adaptation  history  for  this  case.  The  actual  error  is  computed 
relative  to  the  drag  from  a  p  =  3  solution  on  a  8752  element  mesh  generated  by  uniformly 
refining  (in  the  meshing  space)  the  finest  mesh  generated  by  adaptation  on  the  p  =  3 
discretization.  This  value  is  81.344140792  counts.  The  figure  shows  that,  in  general,  the 
p  =  3  discretization  is  slightly  more  efficient,  in  terms  of  DOF  required  to  achieve  a  given 
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(a)  Adapt  Iter  0,  p  —  1,  DOF  =  702 


(b)  Adapt  Iter  2,  p  =  1,  DOF  =  7929 


(c)  Adapt  Iter  0,  p  =  2,  DOF  =  1404 


(d)  Adapt  Iter  2,  p  =  2,  DOF  =  3054 


(e)  Adapt  Iter  0,  p  =  3,  DOF  =  2340 


(f)  Adapt  Iter  2,  p  =  3,  DOF  =  3520 


Figure  5-6:  Skin  friction  distributions  for  flow  over  a  flat  plate  computed  using  the  standard 
weighting  discretization 
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(a)  Adapt  Iter  0,  p  —  1,  DOF  =  702 


(b)  Adapt  Iter  2,  p  =  1,  DOF  =  8850 


(c)  Adapt  Iter  0,  p  =  2,  DOF  =  1404 


(d)  Adapt  Iter  2,  p  =  2,  DOF  =  3390 


(e)  Adapt  Iter  0,  p  =  3,  DOF  =  2340 


(f)  Adapt  Iter  2,  p  =  3,  DOF  =  3610 


Figure  5-7:  Skin  friction  distributions  for  flow  over  a  flat  plate  computed  using  the  dual 
consistent  discretization 
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(a)  Adapt  Iter  0,  p  —  1,  DOF  =  702 


(b)  Adapt  Iter  2,  p  =  1,  DOF  =  8544 


(c)  Adapt  Iter  0,  p  =  2,  DOF  =  1404 


(d)  Adapt  Iter  2,  p  =  2,  DOF  =  3096 


(e)  Adapt  Iter  0,  p  =  3,  DOF  =  2340 


(f)  Adapt  Iter  2,  p  =  3,  DOF  =  3530 


Figure  5-8:  Skin  friction  distributions  for  flow  over  a  flat  plate  computed  using  the  mixed 
formulation  discretization 
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(a)  Meshing  space 


(b)  Physical  space 


Figure  5-9;  Initial  mesh  (1600  elements)  for  the  ellipse  adaptation  test  case,  shown  in  both 
the  meshing  and  physical  spaces 
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(a)  Error  estimate 


(b)  Actual  error 


Figure  5-10:  Estimated  and  actual  drag  error  versus  DOF  for  adaptation  on  the  drag  on  an 
ellipse 


accuracy,  than  the  p  =  2  discretization.  Moreover,  both  p  =  2  and  p  =  3  are  significantly 
more  efficient  that  p  =  1.  For  example,  the  p  =  1  discretization  requires  21999  DOF  to 
achieve  an  error  of  approximately  3.0  x  10“^  counts.  Alternatively,  the  p  =  3  discretization 
gives  an  error  of  approximately  2.1  x  10~^  counts  using  only  6280  DOF. 

Figure  5-11  shows  the  pressure  coefficient  distribution  on  the  initial  and  final  meshes 
for  p  =  1,  2, 3.  While  the  differences  in  the  pressure  coefficient  between  the  initial  and  final 
meshes  are  fairly  minor,  the  skin  friction  changes  dramatically.  Figure  5-12  shows  the  skin 
friction  distribution  on  the  initial  and  final  meshes  for  p  =  1,2,  3.  On  the  initial  mesh,  the 
skin  friction  is  highly  oscillatory  for  all  p.  Furthermore,  the  values  over  much  of  the  airfoil 
are  significantly  larger  on  the  initial  mesh  than  the  final,  leading  to  drag  values  that  are  too 
large.  Despite  the  fact  that  the  initial  solutions  are  severely  under-resolved,  the  adaptation 
algorithm  is  able  to  improve  the  quality  of  the  meshes  and  the  resulting  solutions.  Thus,  on 
the  final  meshes,  while  some  oscillations  are  still  visible,  their  amplitude  has  been  reduced, 
especially  for  p  =  2, 3. 

Finally,  to  demonstrate  the  changes  the  adaptation  scheme  makes  to  the  mesh,  Figure  5- 
13  compares  the  final  mesh  resulting  from  the  p  =  3  run  with  the  initial  mesh.  While  the 
initial  mesh  is  somewhat  anisotropic,  it  is  clear  that,  as  expected,  the  final  adapted  mesh 
has  much  more  resolution  near  the  wall. 
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(a)  j3=  1,  _DOF  =  4800 


(b)  p=l,  DOF  =  21999 


(c)  p  =  2,  DOF  =  9600 


(d)  p  =  2,  DOF  =  31074 


(e)  p  =  3,  DOF  =  16000 


(f)  p  =  3,  DOF  =  21880 


Figure  5-11:  Pressure  coefficient  distributions  on  initial  and  final  meshes  for  flow  over  an 
ellipse 
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(a)  p=l,  DOF  =  4800  (b)  p  =  1,  DOF  =  21999 


(c)  p  =  2,  DOF  =  9600 


(d)  p  =  2,  DOF  =  31074 


(e)  p  =  3,  DOF  =  16000 


(f)  p  =  3,  DOF  =  21880 


Figure  5-12:  Skin  friction  distributions  on  initial  and  final  meshes  for  flow  over  an  ellipse 


119 


5.2.3  NACA  0012 


The  third  test  geometry  is  a  NACA  0012  airfoil.  The  airfoil  definition  and  computational 
domain  for  this  case  are  given  in  Section  3.4.2.  Two  flow  conditions  are  considered:  Moo  = 
0.25,  Rbc  =  1  X  10®,  a  =  0  and  =  0.25,  Rcc  =  1  x  10^,  a  =  0.  In  both  cases,  the  airfoil 
surface  is  a  no  slip,  adiabatic  wall  while  the  outer  boundary  is  treated  with  the  farfield, 
full  state  approach.  Only  the  linear  elasticity  mesh  movement  approach  with  isoparametric 
elements  is  used  in  this  case.  Finally,  the  assumed  error  convergence  rate  was  lowered  from 
s  +  t  =  2p  io  s  +  t  =  p  +  1.  This  step  was  taken  to  improve  the  performance  of  the  scheme 
for  the  initial  adaptation  iterations,  where  the  solution  is  severely  under-resolved.  In  this 
regime,  it  is  unrealistic  to  expect  to  achieve  the  optimal  asymptotic  order  of  accuracy.  Thus, 
the  assumed  rate  was  lowered,  causing  the  adaptation  to  be  more  aggressive — i.e.  decrease 
h  more  rapidly. 

Figure  5-14  shows  the  adaptation  history  for  the  Rec  =  1  x  10®  case.  The  initial  mesh 
for  this  case  was  generated  by  taking  two  p  =  3  adaptation  iterations  on  Moo  =  0.25,  Re  = 
1  X  10®,  a  =  0  flow,  starting  from  a  1280  element  structured  mesh.  The  actual  error  is 
computed  relative  to  the  drag  from  a,  p  =  3  solution  on  a  8170  element  mesh  generated  by 
taking  an  additional  p  =  3  adaptation  iteration  from  the  finest  p  =  3  result  shown.  For  this 
solution,  the  drag  value  is  107.17549806  counts.  As  expected,  the  p  =  2  and  p  =  3  schemes 


(a)  Error  estimate 


(b)  Actual  error 


Figure  5-14:  Estimated  and  actual  drag  error  versus  degrees  DOF  for  adaptation  on  the 
drag  on  a  NACA  0012  in  Rcc  =  1  x  10®  flow 

are  both  more  efficient  than  the  p  =  1  discretization.  For  example,  to  achieve  an  error  of 
7.5  X  10“^  counts,  the  p  =  3  discretization  uses  12110  DOF,  while  the  p  =  1  discretization 
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would  require  more  than  20000. 

The  pressure  coefficient  distribution  on  both  the  initial  and  final  meshes  for  all  p  is 
shown  in  Figure  5-15.  Oscillations  on  the  initial  mesh  are  more  apparent  here  than  in 
the  ellipse  case.  However,  after  adaptation,  smooth  Cp  profiles  are  obtained  for  all  p.  The 
skin  friction  is  also  oscillatory  on  the  initial  mesh,  as  shown  in  Figure  5-16.  Again,  the 
adaptation  algorithm  is  successful  in  improving  the  mesh  and  solution  qualities  despite  the 
poor  quality  of  the  initial  solutions.  Thus,  the  profiles  on  the  final  meshes  are  significantly 
smoother  than  on  the  initial  meshes. 

Figure  5-17  shows  the  adaptation  history  for  the  Rcc  =  1  x  10^  case.  The  initial  mesh  for 
this  case  is  the  mesh  generated  after  four  p  =  3  adaptation  iterations  on  the  Rcc  =  1  x  10® 
case.  The  actual  error  is  computed  relative  to  the  extrapolated  p  =  5  result  given  in 
Section  3.4.2.  For  this  case,  the  p  =  2  method  produces  the  best  results  for  errors  larger 
than  approximately  1  x  10“^  counts.  This  fact  results  at  least  partially  from  the  mesh 
generation  parameters  used  here.  The  BAMG  mesh  generator  bounds  the  rate  of  change  of 
the  mesh  size  by  a  user  specified  parameter  (see  “-ratio”  in  the  BAMG  User’s  Guide  [58]). 
When  “-ratio”  is  set  to  r,  the  rate  of  change  of  the  mesh  spacing  is  bounded  by  log(r).  For 
the  p  =  1,2  cases,  r  =  4  was  used.  For  p  =  3,  the  solution  diverged  after  a  single  adaptation 
iteration  when  using  r  =  4.  Setting  r  =  2  alleviates  this  problem,  but  results  in  larger, 
less  efficient  meshes.  Despite  this  fact,  for  small  error  tolerances,  p  =  3  is  competitive  with 

p  =  2. 

Figure  5-18  shows  the  pressure  coefficient  distribution  on  the  initial  and  final  meshes 
for  p  =  1,2,3.  The  oscillations  observed  on  the  initial  mesh  are  virtually  eliminated  by 
the  adaptive  algorithm.  Figure  5-19  shows  the  skin  friction  distributions.  As  observed  in 
the  previous  cases,  the  oscillations  in  the  skin  friction  are  much  more  severe  than  those  in 
the  pressure  coefficient.  In  the  p  =  2,3  cases,  these  oscillations  are  largely  eliminated  by 
adaptation.  However,  the  DOF  required  to  eliminate  the  oscillations  in  the  skin  friction  is 
larger  than  that  necessary  to  achieve  engineering  accuracy  for  the  drag.  For  p  =  1,  the  skin 
friction  on  the  final  adapted  mesh  is  still  quite  oscillatory.  It  is  expected  that  additional 
adaptation  iterations  would  produce  a  smooth  skin  friction  distribution  if  enough  DOF  were 
allowed. 
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(a)  p=  1,  DOF  =  4428 


(b)  p=  l,DOF  =  28533 


o 

I 


o 

I 


(c)  p  =  2,  DOF  =  8856 


(d)  p  =  2,  DOF  =  75846 


o 

I 


o 

I 


(e)  p  =  8,  DOF  =  14760 


(f)  p  =  3,  DOF  =  40820 


Figure  5-15:  Pressure  coefficient  distributions  on  initial  and  final  meshes  for  Rcc  =  1  x 
flow  over  a  NACA  0012 
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(a)  p=l,  DOF  =  4428  (b)  p  =  1,  DOF  =  28533 


(c)  p  =  2,  DOF  =  8856  (d)  p  =  2,  DOF  =  75846 


(e)  p  =  3,  DOF  =  14760  (f)  p  =  3,  DOF  =  40820 

Figure  5-16:  Skin  friction  distributions  on  initial  and  final  meshes  for  Rec  =  1  x  10®  flow 
over  a  NACA  0012 
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(a)  Error  estimate 


(b)  Actual  error 


Figure  5-17:  Estimated  and  actual  drag  error  versus  DOF  for  adaptation  on  the  drag  on  a 
NACA  0012  in  Rcc  =  1  x  10^  flow 
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(a)  p  =  1,  DOF  =  6438  (b)  p  =  1,  DOF  =  56664 


(c)  p  =  2,  DOF  =  12876  (d)  p  =  2,  DOF  =  152082 


(e)  p  =  3,  DOF  =  21460  (f)  p  =  3,  DOF  =  126170 


Figure  5-18:  Pressure  coefficient  distributions  on  initial  and  final  meshes  for  Rcc  =  1  x 
flow  over  a  NACA  0012 
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(c)  p  =  2,  DOF  =  12876 


(d)  p  =  2,  DOF  =  152082 


x/c 


(e)  p  =  3,  DOF  =  21460 


(f)  p  =  3,DOF  =  126170 


Figure  5-19:  Skin  friction  distributions  on  initial  and  final  meshes  for  Rcc  =  1  x  10^  flow 
over  a  NACA  0012 
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Chapter  6 


Unsteady  Adaptation  Algorithm 


This  chapter  describes  a  modification  to  the  adaptation  algorithm  presented  in  Chapter  5. 
The  new  algorithm,  referred  to  as  the  unsteady  adaptation  algorithm,  is  intended  to  improve 
the  robustness  of  the  standard  procedure  by  removing  the  need  to  converge  a  steady  state 
solution  before  adapting.  The  motivation  for  this  algorithm  is  described  in  more  detail 
in  Section  6.1.  The  algorithm  is  shown  in  Section  6.2,  and  Section  6.3  contains  sample 
numerical  results. 

6.1  Motivation 

One  drawback  of  the  adaptation  algorithm  presented  in  Chapter  5  is  the  need  to  obtain  a 
steady  state  solution  prior  to  adapting  the  mesh.  This  need  places  a  very  high  robustness 
requirement  on  the  flow  solver.  In  particular,  the  solver  must  be  able  to  obtain  the  steady 
flow  and  adjoint  solutions  on  arbitrarily  bad  meshes.  This  requirement  is  particularly 
restrictive  for  high-order  methods,  which  tend  to  produce  oscillations  in  state  variables 
in  under-resolved  regions  of  the  flow.  These  oscillations  often  cause  the  discrete  solution 
to  diverge.  See  [79]  for  an  example  of  this  phenomenon  in  laminar  flow.  Thus,  while 
it  is  important  to  maximize  the  robustness  of  the  solver,  another  method  for  improving 
robustness  is  to  remove  the  requirement  that  a  steady  state  solution  be  obtained  before 
adapting  the  mesh.  One  idea  for  removing  this  requirement  is  to  march  the  solution  forward 
in  time  and  adapt  the  mesh  at  each  time  step. 

6.2  Unsteady  Algorithm 

Underlying  the  unsteady  algorithm  is  a  procedure  for  estimating  the  error  due  to  spatial 
discretization  at  the  current  time  step.  This  error  estimate  is  presented  in  Section  6.2.1. 
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Section  6.2.2  describes  the  unsteady  adaptation  algorithm  for  uniform  time  step.  Finally, 
in  Section  6.2.3,  the  algorithm  is  extended  to  allow  the  time  step  to  vary  throughout  the 
mesh. 


6.2.1  Error  Estimation 


In  the  unsteady  case,  the  error  estimation  algorithm  is  slightly  modified  from  the  standard 
case.  The  output  of  interest  is  taken  to  be  the  average  over  the  current  time  step  of  the 
steady  state  output  of  interest.  Thus,  if  the  output  of  interest  of  the  steady  simulation  is 
denoted  J ,  the  output  of  interest  for  the  ith  time  step  is  given  by 

where  At  =  Only  backward  Euler  time  discretization  is  considered.  Thus,  the 

discrete  solution  from  to  t*  is  constant  in  time — i.e.  u^(x,  t)  =  — which  implies 

that 

At  each  time  step,  the  semi-linear  form  for  the  unsteady  problem  is  given  by 

vjj  {u^H  -  +  Rh{uh,^h), 


Ru,h{uh,vh)  =  ^  X] 


where  Rh  is  the  semi-linear  form  corresponding  to  the  steady  problem.  Thus,  the  discrete 
problem  for  a  single  time  step  takes  the  following  form;  given  G  Vh,  hud  G  Vh 

such  that 

Ru,h{u’‘h,vh)  =  0,  Vvij  G  Vh- 

Moreover,  one  can  view  this  discrete  problem  as  a  spatial  discretization  of  the  following 
temporally  discrete  but  spatially  continuous  problem;  given  G  V,  find  u*  G  V  such 
that 


1 


(u*  -  u*-i)  +  R{u\  v)  =0,  Vv  G  V, 


(6.1) 


Jn 

where  R{-,  ■)  is  the  weak  form  of  the  steady  governing  equations.  Then,  by  the  same  analysis 
as  in  Section  4.1,  the  error  can  be  written  in  terms  of  the  primal  or  the  adjoint  residual; 


Ju{n^)  -  Ju{u^h) 
Ju{n^)  - 
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where  G  Vh  is  arbitrary,  if)  solves  the  adjoint  problem 


V, -ip)  =  -R«,h(u\ uj:^;  V,  tp)  -  JuiuPu'fj;^)  =  0,  Vv  G  +  V, 


and  the  mean-value  linearizations  are 


Ru,h{u\u^h;v,w)  =  X]  /  +  Rh{u\uh-,v,w), 

^  K-&Th 


Of  course,  the  exact  solution,  u*  and  adjoint,  tp,  are  unknown.  However,  making  the 
same  approximations  as  in  the  steady  case,  the  contribution  to  the  error  due  to  spatial 
discretization  of  a  single  element  can  be  estimated  by 

=  2  (\Ru,hi'^H^  i'^h  ~  '^h)\h)\  +  \Rt,h(^H'^  ('^/i  “  h)\J  ■  (6-2) 

Then,  a  global  error  estimate  is  given  by  e  =  Y^kgTh  set  to  the  discrete 

adjoint.  That  is,  tp^j  solves  the  following  problem:  find  'tpu  £  Vh  such  that 

R-u^hI'^hK'^  H , 'tp  h)  ^  I  +  J''[uj^](vH)  =  0,  VvH  G  Vh. 

The  surrogate  solutions  'ipi^  and  Uh  are  constructed  by  taking  element-block  Jacobi  iterations 
on  the  {p  +  l)th  order  unsteady  problem.  This  procedure  is  exactly  analogous  to  that 
presented  in  Section  4.2  for  the  standard  case. 


6.2.2  Comparison  with  Standard  Algorithm 

To  best  understand  the  unsteady  procedure,  it  is  compared  with  the  standard  algorithm. 
Recall  the  form  of  the  standard  adaptation  algorithm,  presented  in  Section  5.1.  The  main 
idea  of  the  unsteady  algorithm  is  to  estimate  the  error  and  adapt  the  mesh  at  every  time  step. 
The  solution  procedure  utilized  in  step  2  of  the  standard  algorithm  reduces  to  Newton’s 
method  as  At  oo.  For  finite  At,  the  state  update  is  that  resulting  from  a  single  Newton 
iteration  on  a  backward  Euler  discretization  of  the  analogous  unsteady  problem.  By  taking 
additional  Newton  iterations  on  this  unsteady  problem,  one  can  find  the  update  that  solves 
the  unsteady  problem,  as  shown  in  Chapter  2.  Given  this  solution  and  the  appropriate 
adjoint,  the  error  at  the  given  time  step  can  be  estimated,  as  demonstrated  in  Section  6.2.1. 
This  error  estimate  motivates  the  following  algorithm: 

1.  Generate  an  initial  mesh  and  initial  condition,  U®.  Initialize  the  CFL  number,  and 
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compute  the  corresponding  time  step,  At.  Set  i  =  1. 

2.  While  ||R(U*-i)|| 

(a)  Solve  the  unsteady  problem  at  the  current  time  step,  Ru(U*)  =  0,  where 

R«(U')  =  -  U'-i)  +  R(U'). 

In  particular, 

i.  Set  WO  =  and  j  =  1. 

ii.  While  ||R„(W-^“^)||  >  6u,toh  where  6u,toi  is  the  unsteady  solution  tolerance, 
and  j  <  Jmax,  where  Jmax  is  the  maximum  number  of  iterations  allowed  to 
solve  the  unsteady  problem, 

A.  Update  the  primal  state:  W-^  =  —  ^R„(W-^“^). 

B.  j  ^  j  +  1. 

hi.  If  ||R„(WJ-^)||  >  ^u,toi,  decrease  the  CFL  number  by  the  user  specified 
factor,  recompute  the  time  step,  and  return  to  step  (i.).  Otherwise,  U*  = 

(b)  Find  'I'  that  solves  the  appropriate  dual  problem  for  the  current  time  step: 

/  1  ^  ^  T  ^ 

— M  H - = - 

5U  5U  u. 

(c)  Use  U*  and  dt  to  compute  the  error  estimate,  e,  for  the  current  time  step. 

(d)  Increase  the  CFL  number  by  the  user  specified  factor,  and  recompute  the  time 
step. 

(e)  If  €min  <  c  <  etoi,  do  not  adapt  the  mesh,  i  i  +  1,  and  return  to  step  (2). 
Otherwise,  adapt  the  mesh,  i  ^  i  +  1,  and  return  to  step  (2). 

3.  If  e  <  etoh  quit.  Otherwise,  the  steady  solution  has  been  obtained,  but  it  is  not  as 
accurate  as  requested.  Thus,  run  standard  adaptation  algorithm. 

The  time  step  in  this  algorithm  is  computed  from  the  CFL  number,  as  shown  in  Section  2.4. 
Furthermore,  the  CFL  number  is  controlled  by  user  specified  increase  and  decrease  factors. 
In  all  cases  shown,  the  CFL  increase  factor  was  set  to  1.5,  and  the  decrease  factor  was  set  to 
2.0.  Finally,  given  the  error  estimate  and  the  surrogate  solution  used  to  set  the  anisotropy, 
the  adaptation  mechanics  used  in  the  unsteady  algorithm  are  the  same  as  those  described 
in  Section  5.1. 
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6.2.3  Extension  to  Variable  Time  Step 


The  algorithm  presented  in  Sections  6.2.1  and  6.2.2  uses  a  time  step  that  is  constant  through¬ 
out  the  mesh.  Thus,  although  only  first-order  accurate,  the  backward-Euler  temporal  dis¬ 
cretization  is  at  least  a  consistent  discretization  of  the  unsteady  problem.  This  algorithm 
should  be  quite  robust — i.e.  one  expects  that  if  small  enough  time  steps  are  used,  the  al¬ 
gorithm  could  successfully  march  through  any  transients  encountered  between  the  initial 
condition  and  steady  state  solution.  However,  given  that  the  time  step  is  set  by  taking 
a  minimum  over  the  entire  mesh,  the  time  step  will  be  restricted  to  that  required  by  the 
smallest  elements.  To  march  to  the  steady  state  solution  more  rapidly,  one  can  consider 
allowing  the  time  step  to  vary  from  element  to  element.  In  this  case,  each  element,  k,  uses 
the  time  step  as  defined  in  Section  2.4.  Thus,  the  semi-linear  form  for  the  fully  discrete 
problem  becomes 


kGTh 


At, 


—  Vh). 


For  this  case,  it  is  not  clear  that  the  corresponding  temporally  discrete  but  spatially  contin¬ 
uous  problem — as  shown  for  uniform  time  step  in  (6.1) — is  well-posed.  However,  in  practice, 
this  is  of  little  concern,  as  one  can  instead  estimate  the  error  relative  to  a  finer  spatial  dis¬ 
cretization.  Specifically,  the  DWR  method  is  applied  to  estimate  the  error  J{vl\)  — 
where  vl\  G  Vh  satisfies 


=  0, 


Vvft  G  Vh, 


and,  for  example,  Vh  is  a  higher-order  space  than  Vh-  The  resulting  error  estimate  has  the 
same  form  as  that  given  in  (6.2).  Moreover,  the  resulting  unsteady  algorithm  is  the  same 
as  that  given  in  Section  6.2.2  except  for  the  use  of  a  different  time  step  on  each  element. 


6.3  Numerical  Results 

To  demonstrate  the  unsteady  algorithm,  it  is  applied  to  three  test  cases.  For  the  first  two 
cases,  simple  geometries  that  have  already  been  considered  in  Chapter  5  are  used.  Unlike 
the  results  of  Chapter  5  however,  the  initial  meshes  used  here  are  obviously  inappropriate 
for  RANS  calculations.  The  final  test  geometry  is  a  multi-element  airfoil.  This  case  is 
started  from  an  inviscid-style,  isotropic  mesh,  and  the  unsteady  algorithm  is  used  to  drive 
the  adaptation  until  the  mesh  is  reasonable  for  obtaining  the  steady  flow  solution. 
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6.3.1  Flat  Plate 


The  first  test  case  is  Rcc  =  1  x  10^,  =  0.25  flow  over  a  flat  plate  with  zero  pressure 

gradient.  The  steady  adaptation  algorithm  is  applied  to  this  case  in  Section  5.2.1.  The 
initial  mesh  used  in  that  case  contained  234  elements,  and  the  height  of  the  first  layer  of 
elements  was  approximately  =  9.5  at  x/c  =  1.0.  The  initial  mesh  used  here,  shown 
in  Figure  6-1,  has  only  78  elements.  Furthermore,  while  the  spacing  in  the  x  direction  is 


Figure  6-1:  Initial  mesh  for  unsteady  adaptation  applied  to  flow  over  a  flat  plate 


the  same  as  the  234  element  mesh,  the  y  spacing  is  much  coarser.  In  particular,  the  height 
of  the  first  layer  of  elements  is  approximately  Ay+  =  2.1  x  10^  at  x/c  =  1.0.  This  spacing 
is  clearly  too  large  for  accurate  RANS  boundary  layer  computations.  More  importantly, 
it  has  not  been  possible  to  obtain  converged  p  >  1,  steady  state  solutions  on  this  mesh. 
However,  the  unsteady  adaptation  algorithm  succeeds  using  this  very  coarse  initial  mesh. 
Specifically,  initializing  the  state  to  uniform  flow,  i.e. 


R-  =  1.0, 

poo 


pu 


Poo^c 


=  1.0, 


pv 


P  00^00 


=  0.0, 


pE 

Pco'u'/x^ 


=  29.07,  —  =  1.0 

Poo 


with  the  initial  CFL  set  to  1.0  and  using  the  p  =  3,  dual  consistent  discretization,  the 
unsteady  adaptation  algorithm  eventually  converges  to  a  steady  state  solution  that  satisfies 
the  desired  error  tolerance  of  0.2  drag  counts. 

The  adaptation  histories  for  both  the  uniform  and  variable  time  step  versions  of  the 
algorithm  are  shown  in  Figure  6-2.  Examining  the  uniform  time  step  results,  the  figure 
shows  that  the  estimated  error  drops  below  the  desired  tolerance  at  iteration  21.  However, 
since  the  spatial  residual  is  not  yet  converged,  the  algorithm  continues  to  march  forward 
in  time  without  adapting  the  mesh.  At  iteration  25,  the  error  estimate  drops  below  the 
minimum  error  request  (one-tenth  of  the  specified  error  tolerance).  Thus,  mesh  adaptation 
is  performed.  In  fact,  the  mesh  is  coarsened  because  the  error  is  smaller  than  required. 
At  iteration  42,  the  error  estimate  re-enters  the  desired  range  for  the  final  time.  Thus,  no 
further  adaptation  is  performed.  Also  at  iteration  42,  the  time  step  is  decreased  by  nearly 
two  orders  of  magnitude.  This  decrease  was  required  to  obtain  a  converged  solution  for  the 
time  step.  The  steady  state  solution  is  obtained  after  an  additional  26  time  steps.  The 
results  for  the  variable  time  step  case  are  qualitatively  similar  except  that  fewer  iterations 
are  required.  In  particular,  the  variable  time  step  algorithm  obtains  the  steady  state  solution 


134 


CFL  Number  Drag  (counts)  DOF 


(a)  DOF  (b)  Error  estimate 


(c)  Drag 


(d)  Spatial  residual  norm 


(e)  CFL  Number  (f)  Time  step  (Uniform  At  only) 


Figure  6-2:  Unsteady  adaptation  history  for  the  flat  plate  test  case 
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in  48  iterations,  compared  to  68  for  the  uniform  time  step  case. 

Meshes  generated  by  the  uniform  time  step  algorithm  are  shown  in  Figure  6-3.  Over  the 
course  of  the  adaptation  run,  the  number  of  DOF  increases  by  a  factor  of  approximately 
2.7  from  the  initial  to  the  final  mesh.  The  largest  mesh  generated  during  the  run  is  larger 
than  the  final  mesh  by  a  factor  of  approximately  1.4. 


(a)  Finest  mesh,  288  elements 


(b)  Final  mesh,  208  elements 


Figure  6-3:  Finest  and  final  meshes  for  unsteady  adaptation  using  uniform  time  step  applied 
to  flow  over  a  flat  plate 

The  final  mesh  generated  using  the  variable  time  step  approach  is  shown  in  Figure  6-4. 
While  this  mesh  has  more  elements  than  the  final  mesh  produced  by  the  uniform  time  step 


Figure  6-4:  Final  mesh  (252  elements)  for  unsteady  adaptation  using  variable  time  step 
applied  to  flow  over  a  flat  plate 

case,  it  also  gives  a  smaller  error  estimate,  although  both  solutions  satisfy  the  requested 
tolerance. 

To  ensure  that  the  unsteady  algorithm  is  producing  reasonable  results,  the  final  error, 
error  estimate,  and  DOF  can  be  compared  to  the  results  of  the  standard  algorithm.  At 
the  converged  steady  state,  the  result  from  the  unsteady  algorithm  with  uniform  time  step 
has  2080  DOF.  The  maximal  error  estimate  (max(ep,.im)  ^duai))  for  the  solution  on  the  final 
mesh  is  2.45  x  10“^  counts,  and  the  actual  error  is  2.94  x  10“^  counts.  Using  the  variable 
time  step  version  of  the  algorithm,  the  final  steady  solution  with  2520  DOF  has  a  maximal 
error  estimate  of  1.37  x  10“^  counts  and  an  actual  error  of  2.31  x  10“^  counts.  Alternatively, 
after  a  single  iteration,  the  standard  algorithm  for  p  =  3,  started  with  2340  DOF,  gives  2630 
DOF  with  a  maximal  error  estimate  of  1.55  x  10“^  and  an  actual  error  of  4.15  x  10“^. 


136 


6.3.2  Ellipse 

The  second  test  case  is  Rbc  =  2  x  10^,  M^o  =  0.25,  a  =  0  flow  over  an  ellipse.  The  geometry 
for  this  case  is  the  same  as  that  considered  in  Section  5.2.2,  but  the  Reynolds  number  is 
twice  as  large.  The  p  =  3,  dual  consistent  discretization  is  used.  The  initial  condition  is 
the  same  as  that  used  for  the  flat  plate  case,  but  the  initial  CFL  was  set  to  0.1.  For  this 
case,  only  uniform  time  step  results  are  presented. 

As  in  the  flat  plate  case,  the  initial  mesh,  shown  in  Figure  6-5,  is  extremely  coarse. 
It  contains  400  q  =  3  elements.  All  meshes  for  this  case  are  generated  using  the  global 
mapping  approach  described  in  Section  5.1.3.  In  the  figure,  each  edge  of  the  true  mesh  is 
approximated  by  linear  interpolation  between  the  q  =  3  nodes  on  the  edge. 

The  adaptation  history  is  shown  in  Figure  6-6.  As  with  the  flat  plate  case,  the  mesh  is 
initially  coarsened.  After  this  initial  coarsening,  elements  are  added  near  the  body  until  the 
error  tolerance  is  satisfied  at  21  iterations.  The  error  estimate  drops  below  the  minimum 
error  request  after  24  iterations,  and  the  mesh  is  coarsened  again.  The  final  mesh  adaptation 
occurs  at  the  37th  iteration,  after  which  26  additional  iterations  are  required  to  converge 
the  steady  solution. 

Figure  6-7  shows  the  finest  mesh,  generated  after  21  iterations,  and  the  final  mesh.  The 
finest  mesh  is  significantly  finer  near  the  ellipse.  This  resolution  is  required  to  resolve  the 
very  thin  boundary  layer  that  forms  after  starting  the  calculation  from  uniform  flow.  As 
expected,  the  final  mesh  shows  significant  refinement  in  the  boundary  layer  relative  to  the 
initial  mesh.  However,  the  wake  is  left  somewhat  under-resolved.  It  is  expected  that  for  a 
smaller  error  tolerance,  refinement  would  occur  in  the  wake  as  well. 


6.3.3  Advanced  Energy  Efficient  Transport  Three-Element  Airfoil 

The  final  test  case  is  Moo  =  0.26,  Rec  =  9  x  10®  (based  on  the  cruise  configuration  chord), 
a  =  8°  flow  around  the  advanced  Energy  Efficient  Transport  (EET)  three-element  airfoil. 
For  this  case,  p  =  3,  isoparametric  elements  are  used,  and  the  curved  elements  are  generated 
via  the  linear  elasticity  node  movement  approach.  The  uniform  time  step  version  of  the 
algorithm  is  used,  and  the  output  of  interest  is  drag. 

The  initial  mesh  used  here,  shown  in  Figure  6-8,  is  an  11063  element,  isotropic  mesh. 
This  mesh  may  be  adequate  for  inviscid  flow,  but  it  is  clearly  not  appropriate  for  a  high 
accuracy  RANS  calculations. 

To  begin  the  calculation,  the  flow  is  initialized  to  zero  flow.  Specifically, 


pu 

Poo^oo 


=  0.0, 


^=0.0,  ^ 

poo^cx)  poo'^oo 


26.416,  =  1.0  X  10"^. 

Poo 
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Drag  (counts)  DOF 


(a)  DOF  (b)  Error  estimate 


(c)  Drag  (d)  Spatial  residual  norm 


(e)  Time  step 


Figure  6-6:  Unsteady  adaptation  history  for  the  ellipse  test  case 
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-100  -50  0  50  100 

(a)  Finest  mesh,  3260  elements 


(c)  Finest  mesh,  zoom 


Figure  6-7:  Finest  and  final  meshes  for  unsteac 
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I 


(b)  Final  mesh,  2538  elements 


(d)  Final  mesh,  zoom 


adaptation  applied  to  flow  over  an  ellipse 


Figure  6-8:  Initial  mesh  for  unsteady  adaptation  applied  to  flow  over  the  EET  3-element 
airfoil 

This  initial  condition — as  opposed  to  the  freestream  condition  used  in  the  flat  plate  and 
ellipse  cases — is  used  to  improve  the  robustness  of  the  procedure.  Using  the  freestream 
flow  initial  condition  here  leads  to  a  mesh  generation  failure  due  to  the  very  small,  high 
aspect  ratio  elements  requested  to  resolve  the  extremely  thin  transient  boundary  layer 
encountered  at  the  beginning  of  the  procedure.  Starting  from  the  zero  flow  condition  relieves 
this  transient  somewhat.  Instead,  a  different  transient  is  encountered  due  to  the  mismatch 
between  the  farfield  boundary  conditions,  which  are  consistent  with  the  freestream  flow, 
and  the  initial  condition. 

For  this  case,  the  unsteady  adaptation  algorithm  is  not  used  to  converge  the  hnal  steady 
state  solution.  Due  to  the  large  cost  of  the  current  algorithm  and  the  length  of  the  transient 
encountered  when  starting  from  zero  flow,  it  is  used  only  as  a  start-up  procedure  to  obtain  a 
reasonable  mesh.  In  particular,  the  unsteady  algorithm  is  used  to  march  the  solution  forward 
in  time  approximately  85  freestream  flow  convective  time  scales.  After  this  time,  the  mesh 
shown  in  Figure  6-9  is  obtained.  Clearly,  this  mesh,  which  contains  10332  elements,  has 
significantly  more  resolution  in  the  boundary  layer  than  the  initial  mesh.  Furthermore,  it  is 
possible  to  obtain  the  steady  state  solution  on  this  mesh,  allowing  the  standard  adaptation 
to  be  applied.  After  a  single  standard  adaptation  iteration,  the  mesh  shown  in  Figure  6-10 
is  generated.  This  mesh  contains  11620  elements  and  is  very  similar  to  the  mesh  generated 


141 


Figure  6-9:  Mesh  obtained  by  unsteady  adaptation  applied  to  flow  over  the  EET  3-element 
airfoil 


Figure  6-10;  Mesh  obtained  by  a  standard  adaptation  iteration  for  flow  over  the  EET 
3-element  airfoil 
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by  the  unsteady  algorithm.  However,  it  is  significantly  finer  near  the  leading  edge  of  the 
slat,  as  illustrated  in  Figure  6-11. 

Figure  6-12  shows  the  solution  obtained  on  the  final  mesh.  The  surface  pressure  co¬ 
efficient  is  shown  in  Figure  6-13.  For  this  solution,  the  lift  coefficient  is  ci  =  3.51,  which 
agrees  well  with  experimental  results  and  other  computations  [71,  5,  103].  Finally,  the  drag 
coefficient  is  =  436  counts,  and  the  error  estimate  is  approximately  3  counts.  This  error 
estimate  gives  a  percent  error  of  approximately  0.7%. 

To  summarize,  while  the  unsteady  procedure  was  not  used  to  converge  the  steady  so¬ 
lution  in  this  case,  it  did  enable  a  reasonable  initial  mesh  to  be  automatically  generated 
starting  from  an  isotropic  mesh  which  was  too  coarse  for  RANS  simulations.  This  reasonable 
initial  mesh  then  allowed  a  steady  solution  to  be  obtained  and  used  to  drive  the  standard 
adaptation  algorithm.  Thus,  the  primary  benefit  of  the  unsteady  procedure  for  this  case 
is  the  automated  generation  of  a  reasonable  initial  mesh.  While  this  is  an  important  step, 
future  work  should  look  to  further  refine  and  improve  the  algorithm. 
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(a)  Indication  of  zoom  region  near  slat  leading  edge 


(b)  Before  standard  adaptation 


(c)  After  standard  adaptation 


Figure  6-11;  Boundary  layer  mesh  on  slat  before  and  after  standard  adaptation  for  EET 
3-element  airfoil 
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(a)  Mach  number 


(b)  Normalized  pressure  {p/ paaUla) 


4000 

2000 

0 


(c)  Normalized  SA  working  variable  (yjvoo) 


Figure  6-12:  Finest  solution  {p  =  3,  11620  elements)  for  the  EFT  3-element  airfoil 
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x/c 


Figure  6-13:  Surface  pressure  coefficient  distribution  for  finest  solution  {p  =  3,  11620  ele¬ 
ments)  for  the  EET  3-element  airfoil 
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Chapter  7 


Conclusions  and  Outlook 


7.1  Summary  and  Conclusions 

This  thesis  presents  research  toward  a  high-order,  adaptive  solution  method  for  the  RANS 
equations  coupled  with  the  SA  turbulence  model.  The  main  contributions  of  this  work  are 
an  analysis  of  the  dual  consistency  of  the  DG  discretization  of  the  SA  model  source  terms, 
an  extension  of  output-based  error  estimation  and  adaptation  to  high-order  discretizations 
of  the  RANS  equations  on  curved,  boundary  conforming  meshes,  and  a  new,  unsteady 
adaptation  framework  to  eliminate  the  need  to  obtain  converged  steady  state  solutions 
prior  to  performing  error  estimation  and  mesh  adaptation. 

The  dual  consistency  analysis  finds  that  the  straightforward  standard  weighting  dis¬ 
cretization  of  gradient  dependent  source  terms  results  in  a  dual  inconsistent  discretization. 
However,  a  dual  consistent  discretization  can  be  constructed  by  adding  a  dual  consistency 
correction  to  the  standard  weighting  scheme.  Furthermore,  discretizations  based  on  the 
mixed  formulation  are  shown  to  be  asymptotically  dual  consistent.  The  impact  of  dual  con¬ 
sistency  on  accuracy  is  illustrated  for  a  scalar  model  problem.  Specifically,  the  results  show 
that  the  dual  consistent  and  asymptotically  dual  consistent  discretizations  achieve  optimal 
accuracy  for  the  primal  solution,  adjoint  solution,  and  a  simple  functional  output.  Alter¬ 
natively,  the  standard  weighting  scheme  produces  suboptimal  results  for  some  p.  For  the 
RANS-SA  system,  the  results  for  a  simple  flat  plate  test  case  demonstrate  that  the  standard 
weighting  discretization  produces  oscillatory  adjoint  solutions  where  the  dual  consistent  and 
mixed  formulation  methods  give  smooth  adjoints.  These  oscillations  appear  most  clearly  in 
the  derivatives  of  the  adjoint  solution.  This  result  agrees  with  the  result  from  the  model 
problem  that  the  adjoint  for  the  standard  weighting  discretization  does  not  converge  to  the 
continuous  adjoint  in  the  broken  norm.  Furthermore,  a  grid  refinement  study  for  the 
NACA  0012  airfoil  shows  that,  when  the  error  is  not  dominated  by  the  boundary  layer  edge 
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singularity,  the  dual  consistent  discretization  is  superior  to  the  standard  weighting  method. 

The  method  of  Dual  Weighted  Residuals  is  applied  to  estimate  the  error  in  functional 
outputs  of  interest.  While  this  method  has  appeared  throughout  the  literature,  to  the  au¬ 
thor’s  knowledge,  this  work  represents  its  first  application  to  high-order  DG  discretizations 
of  the  RANS  equations  on  curved  meshes.  The  implementation  here  relies  on  an  element- 
block  Jacobi  iterative  procedure  to  generate  the  surrogates  for  the  exact  primal  and  adjoint 
solutions  appearing  in  the  error  estimate.  The  results  demonstrate  that  the  error  estimates 
are  generally  accurate  for  both  the  standard  weighting  and  dual  consistent  discretizations. 
However,  especially  for  p  =  2,  the  standard  weighting  error  estimates  can  be  significantly 
smaller  than  the  true  error. 

In  addition  to  the  functional  output  error  estimate,  the  error  estimation  procedure 
provides  an  error  indicator  for  each  element  in  the  mesh.  This  elemental  error  is  used  to 
drive  a  mesh  adaptation  procedure.  Specifically,  the  error  estimate  is  used  to  determine 
the  size  of  elements  desired  in  the  adapted  mesh,  while  the  anisotropy  is  set  by  attempting 
to  achieve  equal  interpolation  error  in  all  directions.  This  procedure  has  been  successfully 
applied  to  multiple  high  Re  test  cases,  showing  that,  given  reasonable  starting  solutions, 
the  algorithm  is  capable  of  improving  the  mesh  and  decreasing  the  output  error.  Moreover, 
it  is  observed  that  the  p  =  2,3  discretizations  are  superior  to  the  p  =  1  in  terms  of  DOF 
required  to  achieve  a  desired  error  tolerance. 

Unfortunately,  the  requirement  that  a  converged  solution  must  be  obtained  before  adapt¬ 
ing  is  quite  restrictive.  In  particular,  it  implies  that  one  must  converge  the  steady  solution 
on  the  initial  mesh.  Of  course,  to  minimize  grid  generation  effort,  it  is  beneficial  to  use 
very  coarse  initial  meshes.  Such  coarse  meshes  can  cause  significant  difficulty  for  high-order 
discretizations,  particularly  for  the  RANS  equations.  Thus,  an  unsteady  adaptation  pro¬ 
cedure,  where  the  need  to  converge  steady  state  solutions  prior  to  adapting  is  eliminated, 
has  been  developed.  In  this  unsteady  procedure,  the  error  due  to  spatial  discretization  is 
estimated  at  every  time  step  and  used  to  drive  the  mesh  adaptation  algorithm.  Numeri¬ 
cal  results  demonstrate  that  this  algorithm  has  the  potential  to  significantly  improve  the 
robustness  of  the  adaptation  procedure  when  starting  from  very  coarse  initial  meshes.  In 
particular,  starting  from  meshes  that  are  clearly  inappropriate  for  RANS  computations, 
accurate  steady  state  p  =  3  solutions  are  obtained. 


7.2  Future  Work 

Many  areas  for  future  work  remain.  A  few  ideas  for  further  research  are  described  here. 
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7.2.1  Turbulence  Model  Improvements 

The  modifications  to  the  SA  model  described  in  Chapter  2  are  designed  to  make  the  model 
equation  stable  for  negative  values  of  the  working  variable.  However,  these  modifications  do 
not  address  the  behavior  that  is  often  the  root  cause  of  the  appearance  of  negative  values: 
the  sharp  transition  in  the  eddy  viscosity  at  the  edge  of  the  boundary  layer.  Neglecting 
the  effect  of  the  viscosity  of  the  fluid,  the  first  derivative  of  the  turbulence  model  working 
variable  is  discontinuous  at  the  boundary  layer  edge.  This  discontinuity  strongly  affects  both 
the  primal  and  adjoint  solutions  at  the  boundary  layer  edge.  This  behavior  is  apparent  in  the 
adjoint  derivative  profiles  (Figures  3-8  and  3-9)  shown  in  Chapter  3.  These  discontinuities  at 
the  boundary  layer  edge,  while  present  in  many  turbulence  models,  are  non-physical  [108]. 
More  importantly,  they  produce  oscillations  in  high-order  solutions  that  can  cause  the 
solution  to  diverge.  Modifications  to  the  model  to  increase  the  smoothness  of  the  eddy 
viscosity  in  this  region  could  lead  to  significant  robustness  improvements  for  high-order 
discretizations  of  the  RANS-SA  system. 

7.2.2  Discretization  Robustness  Improvement 

More  generally,  robustness  for  high-order  RANS  discretizations  is  a  major  issue.  While 
robustness  problems  are  sometimes  tied  to  under-resolution  at  the  boundary  layer  edge,  they 
can  also  be  caused  by  under-resolution  in  other  flow  regions  where  the  turbulence  model  is 
not  at  fault,  e.g.  near  the  leading  edge  of  an  airfoil  where  the  boundary  layer  is  still  laminar. 
One  idea  to  address  such  under-resolved  features  is  the  addition  of  artificial  viscosity.  The 
use  of  artificial  viscosity  has  recently  been  developed  for  shock-capturing  for  high-order 
DG  discretizations  of  the  Euler,  Navier-Stokes,  and  RANS  equations  [84,  10].  In  related 
work,  it  has  been  successfully  applied  to  improve  the  robustness  of  RANS  calculations,  even 
for  flows  without  shocks  [76].  This  technique  should  be  further  investigated  for  high-order 
discretizations  of  the  RANS  equations,  particularly  in  combination  with  adaptation. 

7.2.3  Standard  Adaptation  Algorithm  Modifications 

Two  aspects  of  the  adaptive  procedure  could  be  modified  to  improve  the  performance  of 
the  algorithm.  As  noted  in  Section  5.1.1,  the  desired  anisotropy  is  determined  by  a  high- 
order  analog  of  a  Hessian-based  technique.  Specifically,  approximations  of  the  {p  -|-  l)th 
derivatives  of  the  Mach  number  are  used  to  attempt  to  equidistribute  the  interpolation  error. 
One  problem  with  this  technique  is  that  the  anisotropy  determined  by  the  Mach  number 
may  not  be  appropriate  for  all  solution  and  adjoint  components  over  the  entire  domain. 
For  example,  for  airfoil  flows  with  drag  as  the  output  of  interest,  upstream  of  the  airfoil. 
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the  adjoint  varies  rapidly  in  the  direction  normal  to  the  stagnation  streamline.  Generally, 
the  Mach  number  does  not  capture  this  anisotropy,  potentially  leading  to  poor  resolution 
of  the  adjoint,  which  can  adversely  affect  both  the  error  estimate  and  output  accuracy. 
One  potential  method  for  addressing  this  problem  is  to  compute  the  desired  metric  using 
multiple  quantities  and  find  an  appropriate  intersection  of  the  computed  metrics.  This  idea 
was  explored  by  Castro-Diaz  et  al.  [21]  using  the  primal  state  vector  as  the  quantities  of 
interest.  It  should  be  investigated  using  the  adjoint  state  as  well. 

A  second  potential  deficiency  of  the  algorithm  is  the  reliance  on  an  assumed  order  of 
accuracy.  Even  when  the  asymptotic  order  of  accuracy  is  known,  it  is  unlikely  that  this 
order  will  be  achieved  in  the  initial  iterations  of  an  adaptation  run,  when  the  solution  is  not 
well  resolved.  In  this  case,  the  adaptation  scheme  will  give  less  than  the  desired  reduction  in 
the  error  at  every  iteration,  leading  to  more  iterations  being  required  to  achieve  the  desired 
tolerance.  A  larger  drawback  is  that  the  order  of  accuracy  depends  on  the  regularity  of  the 
exact  solution,  which  may  not  be  known.  Or,  even  if  it  is  known,  the  elements  affected  by 
the  singularities  must  be  determined  from  the  discrete  solution.  Thus,  it  would  be  beneficial 
to  remove  the  order  of  accuracy  from  the  adaptive  algorithm.  One  appealing  approach  is  to 
entirely  avoid  the  metric  specification  step,  instead  driving  the  mesh  adaptation  with  error 
directly,  as  described  by  Park  and  Darmofal  [80]. 

7.2.4  Unsteady  Adaptation  Refinements 

The  unsteady  adaptation  algorithm  presented  here  is  but  a  beginning.  Many  issues  remain 
to  be  resolved.  First,  the  cases  shown  in  Chapter  6  were  started  from  a  uniform  flow  initial 
condition.  While  these  cases  were  successful,  this  initial  condition  causes  a  very  violent 
transient  that  the  solver  must  march  through.  Physically  realistic  initial  conditions  should 
be  explored.  Second,  the  tolerances  used  here  were  set  at  the  levels  desired  for  the  final, 
steady  state  result.  Again,  while  this  setting  worked  in  this  case,  it  is  probably  not  the 
most  efficient  selection.  One  obvious  option  is  to  set  the  error  tolerance  to  some  percentage 
of  the  computed  output.  This  tolerance  may  allow  the  algorithm  to  march  through  the 
initial  transient  using  fewer  DOF,  thus  decreasing  the  total  cost.  Third,  the  stopping 
criterion  should  be  improved.  In  the  cases  shown,  a  significant  fraction  of  the  total  cost 
is  allotted  to  fully  solving  each  time  step  and  estimating  the  error  even  after  the  output 
is  changing  very  little  and  the  adaptation  has  stopped.  This  time  might  be  better  spent 
simply  obtaining  the  steady  state  solution.  Fourth,  only  backward  Euler  time  discretization 
has  been  considered.  The  method  should  be  extended  to  high-order  time  discretizations. 
With  high-order  discretization  in  time,  it  may  be  possible  to  move  through  any  transients 
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more  quickly  and  reliably.  Furthermore,  this  extension  may  be  viewed  as  the  first  step 
toward  a  space-time  adaptation  algorithm  for  truly  unsteady  problems,  where  high-order 
time  discretization  will  be  desirable. 
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Appendix  A 


Derivation  of  the  RANS  Equations 


The  compressible  Navier-Stokes  equations  in  conservation  form  are  given  by 
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where  t  denotes  time,  Xi  is  position,  p  is  density,  Ui  is  velocity,  p  denotes  pressure,  e  is 
the  internal  energy  per  unit  mass,  h  is  the  enthalpy  per  unit  mass,  is  the  viscous  stress 
tensor,  g*  is  the  heat  flux  vector,  and  6ij  is  the  Kronecker  delta.  For  Newtonian  fluids,  the 
viscous  stress  tensor  is  given  by 


Tj  —  2^ 


Iduk  \ 
3dxk  dj  ’ 


where  Sij  =  ^  is  the  strain  rate  tensor,  p  is  the  dynamic  viscosity,  and  the 

second  viscosity  is  assumed  to  be  —  f/^-  Furthermore,  the  heat  flux  vector  is  given  by 
Fourier’s  law; 

dT 

where  T  is  the  temperature  and  k  is  the  thermal  conductivity  of  the  fluid. 

Finally,  it  is  assumed  that  the  fluid  is  an  ideal  gas  with  constant  specific  heats.  Thus, 


p  =  pRT,  e  =  c^T,  h  =  CpT, 

where  R  is  the  ideal  gas  constant,  c„  is  the  specific  heat  at  constant  volume,  and  Cp  is  the 
specific  heat  at  constant  pressure. 
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While  (A.l)  through  (A. 3)  govern  the  motion  of  fluids  in  both  laminar  and  turbulent 
regimes,  the  large  range  of  spatial  and  temporal  scales  present  in  turbulent  flows  make 
direct  solution  of  the  compressible  Navier-Stokes  equations  prohibitively  expensive  for  the 
majority  of  flows  of  engineering  interest  [37,  74],  Thus,  it  is  common  to  instead  solve  the 
RANS  equations,  which  govern  the  turbulent  mean  flow. 

Following  Reynolds  [91],  the  flow  variables  are  written  in  terms  of  mean  and  fluctuating 
parts.  Then,  equations  governing  the  mean  flow  are  derived  by  averaging  the  Navier-Stokes 
equations.  A  summary  of  the  relevant  definitions  and  resulting  equations  is  given  here.  For 
more  details,  see  [108]. 

Given  a  flow  variable,  u(x,t),  define  the  time  average,  u(x),  by 

rto+T 

u(x)  =  lim  /  v{'x.,t)dt. 

Then,  the  flow  variables  can  be  decomposed  into  mean  and  fluctuating  parts.  For  example, 
the  density  can  be  written 

P  =  P  +  p', 

where  the  fluctuating  part,  p' ,  is  defined  by  this  relationship. 

In  addition,  for  compressible  flows,  it  is  convenient  to  define  a  density-weighted  time 
average.  The  density- weighted  average  of  a  quantity  u(x,t)  is  denoted  by  u(x),  where 

2  rto+T 

u(x)  =  -  lim  /  dt. 

P  T--00 

This  averaging  procedure  is  known  as  Favre  averaging.  As  with  the  conventional  average, 
the  density-weighted  average  allows  the  flow  variables  to  be  written  in  terms  of  mean  and 
fluctuating  parts.  For  example, 

Ui  =  Ui  +  u'l- 


Applying  the  time  averaging  procedure  to  the  compressible  Navier-Stokes  equations,  one 
obtains  the  following  equations  governing  the  mean  flow: 
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Equations  (A. 4)  through  (A. 6)  contain  mnltiple  correlation  terms — e.g.  pu'-u'l  in  the  mo¬ 
mentum  equation — that  do  not  appear  in  (A.l)  throngh  (A. 3)  and  are  not  known  in  terms 
of  the  mean  flow  variables.  Thus,  there  are  more  unknowns  than  equations,  which  implies 
that  the  system  is  not  closed.  To  close  the  system,  it  is  necessary  to  relate  these  correlation 
terms  to  the  mean  flow  via  closure  approximations. 

In  general,  one  mnst  provide  closnre  approximations  for  five  terms:  the  Reynolds  stress 
tensor,  pu'-u'l]  the  tnrbnlent  kinetic  energy  per  unit  volnme,  pk  =  ^pu'lu'l]  the  turbulent 
heat  flux  vector,  pujh";  the  molecular  diffusion  term,  Tjiu";  and  the  turbulent  transport 
term,  pu'-^u"u”.  In  this  work,  the  turbulent  kinetic  energy,  molecular  diffusion,  and  turbu¬ 
lent  transport  terms  are  ignored.  This  approximation  is  reasonable  whenever  /c  <C  /i,  which 
holds  for  most  engineering  flows  into  the  supersonic  regime  [108].  Thus,  only  two  closnre 
approximations  are  reqnired.  For  the  Reynolds  stress  tensor,  the  Bonssinesq  approximation 
is  applied.  The  Boussinesq  approximation  relates  the  Reynolds  stress  to  the  mean  flow 
viscous  stress  tensor: 


puju'l  —  2pt  ^Sji  2>dxk^^^ 

where  pt  is  the  eddy  viscosity.  The  addition  of  —2pk6ji/3  is  made  to  guarantee  that  the 
trace  of  the  Reynolds  stress  is  —2pk.  However,  in  keeping  with  the  assumption  that  the 
tnrbnlent  kinetic  energy  is  small  compared  to  the  mean  flow  enthalpy,  this  term  is  neglected 
here. 

Then,  following  the  Reynolds  analogy,  the  turbulent  heat  flnx  vector  is  given  by 
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where  Prt  denotes  the  tnrbnlent  Prandtl  nnmber,  which  is  taken  to  be  Prt  =  0.9. 


(A.8) 
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Appendix  B 


Laminar  NACA  0012  Results 


The  test  case  is  M^o  =  0.25,  Rcc  =  1  x  10^,  a  =  0  flow  over  a  NACA  0012  airfoil.  The 
computational  domain  and  meshes  used  for  this  case  are  described  in  detail  in  Section  3.4.2. 

Figure  B-1  shows  the  drag  error  results.  The  error  is  computed  relative  to  an  “exact” 
drag  value,  which  is  computed  by  extrapolating  the  p  =  5  drag  results,  assuming  that 
the  error  convergence  rate  is  the  same  for  both  refinements.  The  resulting  drag  value  is 
375.2615065  counts. 

The  figure  shows  that  p  =  1,2,3  achieve  nearly  optimal  order  of  accuracy — i.e.  — for 
the  second  mesh  refinement.  While  p  =  4, 5  are  not  achieving  their  optimal  rates,  there  is  no 
evidence  that  the  accuracy  will  asymptote  to  Oih^)  as  in  the  RANS-SA  case.  Furthermore, 
by  examining  the  elemental  error  estimates,  it  is  possible  to  confirm  that  the  boundary 
layer  edge  region  is  not  the  source  of  the  suboptimal  p  =  4,  5  results.  Figure  B-2  shows 
the  near-airfoil  p  =  4  elemental  error  estimates  on  the  coarse  and  fine  meshes.  The  figure 
shows  that  on  the  coarse  mesh,  the  error  is  dominated  by  a  few  elements  near  the  edge  of 
the  boundary  layer  and  wake.  However,  as  expected  for  smooth  solutions,  these  errors  are 
effectively  decreased  by  h  refinement  with  the  high-order  discretization,  and  these  regions 
to  not  appear  to  contribute  significantly  to  the  error  on  the  fine  mesh. 

The  p  =  4  elemental  error  estimates  for  entire  domain  are  shown  in  Figure  B-3.  As 
expected,  the  coarse  mesh  errors  in  the  outer  regions  of  the  domain  are  small  compared  to 
those  near  the  airfoil.  However,  for  the  fine  mesh,  the  error  is  largest  on  the  outer  boundary 
at  the  top  and  bottom  points  of  the  circle.  These  points  are  where  the  boundary  condi¬ 
tion  switches  from  inflow  to  outflow  conditions.  Thus,  it  appears  that  it  is  the  boundary 
conditions  which  cause  the  loss  of  accuracy  for  p  =  4,  5  in  this  case. 

To  summarize,  given  that  this  case  is  run  using  the  same  meshes  as  the  turbulent  NACA 
0012  case,  it  reinforces  the  argument  that  the  RANS-SA  solution  behavior  at  the  boundary 
layer  edge  is  responsible  for  the  suboptimal  accuracy  observed  in  Section  3.4.2. 
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(a)  Drag  error  versus  h  for  p  =  1  through  p  =  5  (p  =  1  :  p  =  2  :  o, 

p  =  3;*,p  =  4:<l,  p  =  5:  x) 


(b)  Drag  error  versus  DOF 


Figure  B-1:  Drag  error  for  laminar  flow  over  a  NACA  0012 
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(a)  Coarse  mesh 
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(b)  Fine  mesh 


Figure  B-2;  Elemental  error  estimate  for  the  p  =  A  discretization  on  the  coarse  and  fine 
meshes  for  the  laminar  NACA  0012  test  case  in  the  near  field 


(a)  Coarse  mesh 


(b)  Fine  mesh 


Figure  B-3:  Elemental  error  estimate  for  the  p  =  A  discretization  on  the  coarse  and  fine 
meshes  for  the  laminar  NACA  0012  test  case  for  the  entire  domain 
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Appendix  C 


A  Priori  Output  Error  Estimation 


This  appendix  provides  a  priori  error  bounds  for  functional  outputs  computed  from  fi¬ 
nite  element  solutions.  The  analysis  demonstrates  the  importance  of  dual  consistency  for 
accurately  computing  functional  outputs.  Furthermore,  the  bounds  obtained  are  used  to 
motivate  the  convergence  rates  used  by  the  mesh  optimization  algorithm  discussed  in  Sec¬ 
tion  5.1.2.  For  simplicity,  a  linear  problem  and  linear  functional  output  are  considered  first. 
Then,  the  analysis  is  extended  to  the  nonlinear  case. 


C.l  Linear  Analysis 

Consider  the  following  continuous  primal  problem;  compute  J (u)  where  u  gV  satisfies 

B{u,v)  =  i{v),  Vu  G  V. 

Here,  B{-,  •)  is  a  bilinear  functional  and  i,  and  J  are  linear  functionals.  The  corresponding 
dual  problem  is  given  by  the  following;  find  V’  £  V  such  that 

B{v,'ip)  =  J{v),  VuGV. 

Assume  that  both  the  continuous  primal  and  dual  problems  are  well-posed. 

Furthermore,  consider  the  following  discretization  of  the  primal  problem;  compute 
J{uh)  where  Uh  G  Vh  satisfies 

Bh{uh,Vh)  =  ^{vh),  Vuft  G  Vh-, 

where  Vh  is  an  appropriate  finite  dimensional  space.  Then,  the  discrete  dual  problem  is  as 
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follows;  find  V’/i  £  V/i  such  that 


Vu/i  G  V/j,. 

As  with  the  continuous  problems,  assume  that  both  the  discrete  primal  and  dual  problems 
are  well-posed. 

Then,  using  the  definitions  of  the  continuous  and  discrete  dual  problems,  the  error  in 
the  computed  functional  output  can  be  written  as 

J{u)  -J {uh)  =  B{u,  V’)  -  Bh{uh,  i^h)- 

To  continue,  assume  that,  for  the  exact  solutions  u  and  V’j  B{u,'4))  =  Bh{u,il!).  Then, 

J{u)-J{uh)  =  Bh{u,^/j)  -  Bh{uh,'il^h) 

=  Bh{u,  -  Bh{uh,  'ip)  +  Bh{uh,  Ip)  -  Bh{uh,  iph) 

=  Bh{u-Uh,ip)  +  B{uh,p^-^ph) 

=  Bh{u  -  Uh,  Ip)  -  Bh{u  -  Uh,  iph)  +  Bh{u  -  Uh,  iph)  +  B{uh,  ip  -  ^ph) 

=  Bh{u  -Uh,^^  -  'iph)  +  Bh{u  -  Uh,  'iph)  +  Bh{uh,  Ip  -  iph)- 

Thus,  if  the  discretization  is  consistent  and  dual  consistent, 

J{u)-  J {uh)  =  Bh{u  -Uh,'ip-  iph)- 

Finally,  if  the  bilinear  form  Bh  is  continuous, 

\J{u)  -J{uh)\  <  C\\\u  -  U;,|||  IIIV’  -  iph\\\, 

where  |||  •  |||  is  an  appropriate  energy  norm.  Assuming  that  the  discretization  achieves 
Ilk  -  Uh\\\  <  Ch^  and  \\\ip  -  iph\\\  <  Ch\ 

\J{u)-J{uh)\<Ch^^K 

For  DG  discretizations  of  hyperbolic  problems,  assuming  the  exact  solutions  u  and  ip  are 
sufficiently  regular,  r  =  s  =  p  +  1/2  (see  Johnson  and  Pitkaranta  [63]).  Thus, 

\J{u)-J{uh)\  <Ch^P+\ 

For  elliptic  problems,  assuming  the  exact  solution  is  sufficiently  regular,  the  BR2  scheme. 


162 


as  well  as  many  other  DG  methods,  achieve  r  =  s  =  p  (see  Arnold  et  al.  [7]).  Thus, 

\J{u)-J{uh)\  <Ch^P. 

C.2  Nonlinear  Analysis 

Consider  the  following  continuous  primal  problem;  compute  77 (u)  where  u  gV  satisfies 

R{u,v)  =  i{v),  Vu  G  V, 

where  R{-,  ■)  is  a  semi-linear  (linear  in  the  second  argument)  functional,  7  is  a  linear  func¬ 
tional,  and  77  is  a  nonlinear  functional.  The  corresponding  dual  problem  is  given  by  the 
following:  find  V'  G  V  such  that 

R'[u\{v,i^)  =  J'[u\{v),  VuGV. 

Assume  that  both  the  continuous  primal  and  dual  problems  are  well-posed. 

Furthermore,  consider  the  following  discretization  of  the  primal  problem:  compute 
J{uh)  where  Uh  G  Vh  satishes 


Rh{uh,Vh)  =  7(u/j),  Mvh  G  Vh- 

Then,  define  the  discrete  dual  problem  is  as  follows:  find  V’/i  £  Vh  such  that 

R'h[u]{vhAh)  =  J'[u]{vh),  VvhGVh- 

As  with  the  continuous  problems,  assume  that  both  the  discrete  primal  and  dual  problems 
are  well-posed. 

The  central  idea  of  the  nonlinear  analysis  shown  here  is  the  examination  of  Taylor  series 
expansions  of  the  output  and  residual  forms  about  the  exact  solution,  u.  In  particular, 
consider  the  Taylor  series  expansion  of  the  function  output; 

J{uh)  =  J{u)+  J'[u]{uh  -u)  +  rout{u,Uh), 

where  ^ 

rout{u,Uh)=  /  J"[u  +  6{uh-u)]{uh-u,uh-u){l-6)d6. 

Jo 

For  the  conditions  required  for  this  statement  to  be  valid,  see  [47]. 
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Assuming  that  the  bilinear  functional  J''[u+9{uh—u)\{- ,  ■)  is  continuous  for  all  9  G  (0, 1), 

\rout{u,Uh)\  <  f  \J”[u  +  9{uh  -  u)]{uh  -u,Uh-  u){l  -  9)\ d9  <  C\\\uh  -  u\\\‘^. 

Jo 

Thus, 

\J{uh)  -  J{u)\  <  \  j'[u]{uh  -  w)!  +  C\\\uh  -  u\f. 

To  complete  the  analysis,  it  is  necessary  to  analyze  the  linear  term,  —  u).  By 

the  definitions  of  the  dual  problems, 

J'[u\{uh  -u)  =  Rh[u\{uh,'iph)  -  R'[u]{u,'4!). 

Assuming  that,  for  the  exact  solutions  u  and  V’j  =  R'f^[u\{u,^),  gives 

J'[u]{u-uh)  =  Rh[u]{u,i;)  -  R',^[u]{uh,i^h) 

=  R'f^[u]{u,i^)  -  R',^[u]{uh,i^)  +  R'h[u]{uh,i^)  -  Rh[u]{uh,i^h) 

=  R'h[u]{u +  R'h[u]{uh,i^ -i^h) 

=  Rh[u]{u  -  Uh,^^)  -  R'h[u]{u  -  Uh,i^h)  +  R'hW]{u  -  Uh,^h)  +  Rh[>A{uh,'il^  -  i^h) 
=  R'hbAiu  -uh,i)  -  iph)  +  R'hbAiu  -  uh,'iph)  +  Rh[u]{uh,ip  -  iph)- 

If  the  scheme  is  dual  consistent,  R'f^[u]{uh,ip  —  iph)  =  0.  The  term  R'f^[u\{u  —  Uh,iph)  corre¬ 
sponds  to  the  consistency  error  term  from  the  linear  analysis.  However,  unlike  the  linear 
case,  this  term  is  not  zero,  even  for  a  consistent  discretization.  To  bound  this  term,  begin 
by  writing  a  Taylor  series  for  the  functional  Rh{-,iph)  about  u.  That  is, 

Rh{Uh,'>Ph)  =  Rhiu^lph)  +RhW]{Uh  -U,i^h)  +rres{u,Uh), 


where 


rresiu,Uh)=  /  Rh[u  +  9{uh-u)]{uh-u,uh-u,'tph){'i--9)d9. 
Jo 


Thus,  assuming  that  tri-linear  functional  R'/^[u  +  9{uh  —  •,  •)  is  continuous  for  all  9  G 

(0,1), 

rres{u,Uh)  <  /  \Rh[u  +  9{uh-u)]{uh-u,Uh-u,'iph)\d9  <C\\\uh-u\\ 

Jo 

If  the  discretization  is  consistent,  Rh{uh,'iph)  —  Rh{u,'iph)  =  0.  Thus, 

\R'hW]{uh  -  u,i^h)\  <  C\\\uh  -  will 
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Further  assuming  that  the  bilinear  form  •)  is  continuous, 


\J'[u]{u  -  uu)\  <  C\\\uu  -  u|||  dllV’,  -  ^lll  +  |||n,  -  u|||) . 


Thus, 


\J{uh)  -  J{u)\  <  C\\\uh  -  u|||  (IllV’ft  -  '0111  +  \\\uh 
Finally,  if  the  scheme  achieves  \\\uh  —  u|||  <  Ch^  and  ,\\\'4>h  —  '0111  < 


J'[u]{u-uh)  <  C7/i“W^+h2^). 
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Appendix  D 


Error  Convergence  Rates  for 
Adaptation 


The  results  shown  in  Chapter  5  demonstrate  that  the  p  =  2,3  adaptive  discretizations  are 
superior  to  the  p  =  1  in  terms  of  DOF  required  to  achieve  a  given  accuracy.  However, 
the  results  should  also  be  compared  to  the  optimal  adaptive  performance.  The  optimal 
performance  of  the  adaptive  scheme  is  not  immediately  obvious.  Thus,  to  get  an  idea  of 
what  to  expect,  two  different  regimes  are  considered. 

In  the  first  regime,  all  flow  features  are  relatively  well  resolved.  Specifically,  assume 
that  the  error  is  approximately  equidistributed  throughout  the  domain  and  that  the  mesh 
anisotropy  is  appropriate — i.e.  there  are  no  regions  of  the  mesh  where  decreasing  the  mesh 
spacing  in  one  direction  is  significantly  more  effective  at  driving  down  the  error  than  decreas¬ 
ing  the  spacing  in  an  orthogonal  direction.  Then,  adaptation  and  global  uniform  refinement 
should  provide  the  same  error  reduction  with  DOF.  Thus,  assuming  that  the  optimal  order 
of  accuracy  for  the  output  of  interest  is  0{h'’),  and  noting  that,  for  isotropic  refinement, 
Ndof  oc  h~^,  where  n  is  the  dimension  of  the  domain,  the  expected  convergence  rate  is  given 
by 

7—1  AT - V  /  71 

■ 

In  the  second  regime,  assume  that  the  output  error  is  dominated  by  error  in  a  particular 
region  of  the  mesh.  In  this  case,  one  would  expect  the  error  reduction  with  DOF  to  be 
dictated  by  the  details  of  the  region  where  the  error  is  large.  For  example,  if  the  total  error 
and  DOF  in  the  mesh  are  dominated  by  some  anisotropic  flow  feature — e.g.  a  shock  or  a 
turbulent  boundary  layer — it  is  likely  that  the  grid  spacing  normal  to  this  feature  is  too 
large.  Thus,  one  might  expect  that  the  error  would  achieve  the  optimal  rate  with  respect 
to  the  grid  spacing  normal  to  the  feature.  That  is,  if  hi  is  the  grid  spacing  normal  to  the 
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feature,  E  on  h\.  Then,  if  the  mesh  spacing  is  decreased  only  in  this  direction,  Ndoj  oc 
Thus, 

E  oc  AT- . 

For  shock-free,  attached  flow,  RANS  calculations  on  coarse  meshes — in  particular  meshes 
where  the  near  wall  spacing  is  too  large — one  expects  the  drag  error  to  be  dominated  by 
errors  in  the  boundary  layer  and,  more  specifically,  errors  in  the  laminar  sublayer.  Fur¬ 
thermore,  the  solution  in  the  laminar  sublayer  is  smooth  and  anisotropic.  Thus,  one  might 
expect  E  oc  However,  given  that  E  oc  is  an  asymptotic  convergence  rate,  it  is  not 

clear  that  this  rate  will  be  achieved  until  the  boundary  layer  is  at  least  somewhat  well  re¬ 
solved.  This  fact  may  lead  to  suboptimal  convergence.  On  the  other  hand,  if  other  regions  of 
the  flow  are  over-resolved,  the  adaptation  algorithm  will  coarsen  the  mesh  in  these  regions. 
This  subtraction  of  DOF  in  other  regions  may  lead  to  better  than  expected  performance  in 
terms  of  error  versus  total  DOF. 

Once  the  near  wall  solution  is  well  resolved,  the  error  may  be  dominated  by  errors  in 
the  boundary  layer  edge  region.  Given  that,  for  the  grid  resolution  typical  in  engineering 
calculations,  the  solution  in  this  region  is  singular,  optimal  convergence  in  this  regime  should 
not  be  expected. 

To  examine  the  error  reduction  versus  DOF  obtained  by  the  adaptation  algorithm, 
results  for  the  ellipse  and  NACA  0012  test  cases  are  shown  in  Tables  D.l  through  D.3. 
In  the  tables,  the  column  Rate  refers  to  the  convergence  rate  of  the  error  with  DOF — i.e. 
E  oc  The  notation  Rate  =  NA  indicates  that  either  the  total  DOF  decreased  or  the 

error  increased  at  that  iteration.  A  star  (*)  in  the  column  labeled  “BL  Edge”  indicates  that 
the  error  estimate  in  the  boundary  layer  edge  region  contributes  significantly  to  the  total 
error  estimate.  Specifically,  if  there  is  an  element  in  the  boundary  layer  edge  region  where 
the  elemental  error  estimate  is  equal  to  or  greater  than  one  half  of  the  maximum  elemental 
error  estimate  in  the  mesh,  then  (*)  appears  in  the  BL  Edge  column.  This  data  is  included 
to  give  a  rough  feel  for  whether  the  boundary  layer  edge  is  playing  an  important  role  in  the 
total  error. 

Eor  p  =  1,  the  adaptation  generally  gives  better  than  the  rate  expected  for  isotropic 
refinement,  and  the  boundary  layer  edge  region  is  never  a  significant  contributor  to  the 
total  error.  Eor  p  =  2, 3  it  is  difficult  to  draw  any  firm  conclusions.  In  general,  it  appears 
that  the  initial  adaptation  iterations  perform  well  in  that  they  are  able  to  decrease  the  error 
or  hold  it  fixed  while  decreasing  the  total  DOE.  In  this  regime,  the  elemental  error  estimates 
are  largest  at  the  wall,  and  DOE  are  added  near  the  wall.  The  decrease  in  total  DOE  is 
due  to  coarsening  of  other  regions  of  the  domain  that  initially  contain  more  resolution 
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than  necessary.  As  the  adaptation  progresses,  the  boundary  layer  edge  begins  to  play  a 
significant  role.  However,  it  is  not  clear  that  this  fact  is  responsible  for  the  degradation  of 
the  convergence  rate  observed  in  both  the  ellipse  and  NACA  0012,  Rec  =  1  x  10^^  cases. 
Finally,  it  should  be  pointed  out  that  the  most  consistent  result  in  terms  of  the  observed 
convergence  rate  is  that  for  the  p  =  3,  NACA  0012  at  Rcc  =  1  x  10^,  where  the  BAMG  grid 
spacing  ratio  was  decreased  to  provide  a  smoother  grid. 


Table  D.l;  Error  convergence  versus  DOF  for  the  ellipse  test  case 


p 

Adapt  Iter 

Ndof 

Error  (counts) 

Rate 

BL  Edge 

1 

0 

4800 

3.16el 

— 

1 

1 

13605 

8.69e0 

-1.238 

1 

2 

21999 

3.03e-2 

-3.501 

2 

0 

9600 

2.478el 

— 

2 

1 

6078 

9.870e0 

NA 

2 

2 

5928 

2.871e-2 

NA 

2 

3 

7926 

1.828e-l 

NA 

2 

4 

14904 

6.140e-3 

-5.374 

2 

5 

31074 

1.135e-3 

-2.297 

3 

0 

16000 

1.315el 

— 

3 

1 

8400 

5.164e0 

NA 

3 

2 

4450 

7.686e0 

NA 

3 

3 

4740 

1.755e-l 

NA 

3 

4 

5220 

5.583e-l 

NA 

3 

5 

6280 

2.133e-2 

-17.659 

3 

6 

9470 

3.160e-3 

-4.649 

3 

7 

14310 

2.766e-3 

-0.323 

3 

8 

21880 

1.122e-3 

-2.126 

Given  that  the  results  do  not  agree  with  the  expected  optimal  rates  and  do  not  support 
any  obvious  conclusion,  further  work  is  required  to  more  fully  understand  the  performance 
of  the  algorithm  in  this  respect. 
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Table  D.2:  Error  convergence  versus  DOF  for  the  NACA  0012  {Rec  =  1  X  10®)  test  case 


p 

Adapt  Iter 

^dof 

Error  (counts) 

T 

BL  Edge 

1 

0 

4428 

8.734el 

— 

1 

1 

8115 

7.467e0 

-4.060 

1 

2 

28533 

4.113e-2 

-4.137 

2 

0 

8856 

1.964el 

— 

2 

1 

5670 

1.171e0 

NA 

2 

2 

10536 

3.064e-l 

-2.164 

2 

3 

27168 

9.569e-3 

-3.660 

* 

2 

4 

75846 

8.980e-4 

-2.305 

* 

3 

0 

14760 

6.375e0 

— 

3 

1 

6740 

1.689e0 

NA 

3 

2 

8710 

1.176e0 

-1.411 

3 

3 

12110 

7.501e-2 

-8.351 

3 

4 

21460 

8.155e-4 

-7.903 

* 

3 

5 

40820 

7.803e-4 

-0.069 

* 

Table  D.3;  Error  convergence  versus  DOF  for  the  NACA  0012  {Rcc  =  1  x  10^)  test  case 


P 

Adapt  Iter 

Ndof 

Error  (counts) 

r 

BL  Edge 

1 

0 

6438 

5.872el 

— 

1 

1 

19659 

8.288e0 

-1.754 

1 

2 

56664 

2.430e-l 

-3.334 

2 

0 

12876 

3.054el 

— 

2 

1 

12774 

1.512e0 

NA 

2 

2 

20562 

1.948e-l 

-4.304 

2 

3 

55884 

2.874e-2 

-1.914 

2 

4 

152082 

1.059e-3 

-3.298 

3 

0 

21460 

5.437e0 

— 

3 

1 

30530 

1.471e0 

-3.708 

3 

2 

46490 

1.661e-l 

-5.187 

3 

3 

76380 

1.748e-2 

-4.535 

3 

4 

126170 

1.676e-3 

-4.672 
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Appendix  E 


Boundary  Conditions 


This  appendix  details  the  boundary  conditions  used  in  this  work.  All  boundary  conditions 
are  enforced  weakly  through  the  fluxes,  and  across  the  domain  boundary.  In  general, 
these  fluxes  are  computed  by  evaluating  the  flux  using  boundary  state  and  state  gradient 
vectors  that  depend  on  both  the  interior  state  and  gradient  and  the  boundary  condition 
information.  The  calculation  of  the  boundary  state  and  state  gradient  as  well  as  departures 
from  this  framework  are  described  for  various  boundary  condition  options  below. 

E.l  Farfield,  Full  State  Boundary 

At  a  farfield  “full  state”  boundary,  the  boundary  state  vector,  u^,  is  specified  by  the  user. 
Then,  the  inviscid  flux,  ^  is  computed  using  the  Roe  flux; 

=  ?^(u+,u^n+). 

The  viscous  flux  is  calculated  using  the  user  specified  boundary  state  and  the  state  gradient 
evaluated  from  the  interior: 

Tl  =  A(vt)VMl. 

As  shown  by  Lu  [72],  the  specification  of  the  inviscid  flux  in  this  manner  renders  the 
discretization  dual  inconsistent.  Thus,  this  boundary  condition  is  only  used  for  farfield 
boundaries  where  this  dual  inconsistency  should  have  little  impact  on  the  flow  solution. 

E.2  Subsonic  Inflow:  pt,  a,  pu 

The  subsonic  inflow  conditions  are  set  using  Riemann  invariants.  For  2-D  flows,  at  subsonic 
inflow  boundaries,  analysis  of  the  convective  terms  of  the  RANS-SA  system  shows  that  there 
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are  four  incoming  and  one  outgoing  Riemann  invariants.  Thus,  for  convection  dominated 
flow,  to  appropriately  specify  a  subsonic  inflow  condition,  four  quantities  must  be  specified. 
In  this  work,  the  total  temperature,  Tt,  total  pressure,  pt,  inflow  angle,  a,  and  turbulence 
working  variable  pP  are  specified.  This  information  is  combined  with  the  outgoing  Riemann 
invariant  to  compute  the  boundary  state.  Then,  the  inviscid  flux,  is  given  by 

=  .F(u^). 

The  viscous  flux  is  calculated  using  boundary  state  and  the  state  gradient  evaluated  from 
the  interior, 

=  ^(u')Vu+. 

E.3  Subsonic  Outflow:  p 

The  subsonic  outflow  boundary  condition  is  also  specified  using  Riemann  invariants.  In  this 
case,  there  is  one  incoming  invariant  and  four  outgoing.  Thus,  to  appropriately  specify  a 
subsonic  outflow  condition,  one  quantity  must  be  specified.  In  this  work,  the  static  pressure, 
p,  is  specified.  This  quantity  is  combined  with  the  outgoing  Riemann  invariants  to  compute 
the  boundary  state.  Then,  the  inviscid  flux,  is  given  by 

=  jr(u^). 

The  viscous  flux  is  calculated  using  boundary  state  and  the  state  gradient  evaluated  from 
the  interior, 

=  ^(u')Vu+. 

E.4  No  Slip,  Adiabatic  Wall 

At  a  no  slip,  adiabatic  wall,  the  desired  conditions  are  zero  velocity,  heat  flux,  and  eddy 
viscosity  at  the  wall.  Using  the  this  information,  the  boundary  state  is  given  by 

=  [p+,  0,  0,  pE+,  0]'^, 

and  the  inviscid  flux,  is  given  by 


=  E(u^). 

When  combined  with  the  no  slip  condition,  the  adiabatic  condition  requires  that  the 
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viscous  flux  in  the  energy  equation  is  zero.  Thus,  the  viscous  flux  is  computed  using 
the  interior  gradient,  except  for  the  energy  equation  component,  which  is  set  to  zero.  In 
particular, 


{^v  -  h/i/(u/i))  •  n+  = 


(^^(u'’)Vu+  -  ■  n+ 

(^(u'’)Vu+  -  ■  n+ 

^^(u^)Vu+  -  •  n+ 


0 


E.5  Symmetry  Plane 


The  symmetry  plane  boundary  condition  requires  that  the  state  be  symmetric  about  the 
plane  and  that  its  derivatives  normal  to  the  plane  are  continuous.  These  conditions  require 
that  the  flow  is  tangent  to  the  plane.  Thus,  the  boundary  state  is  set  to  the  interior  state 
with  the  normal  velocity  subtracted: 


=  P'^, 

pu\ 

=  put 

pA 

=  put 

pE^ 

=  pE+ 

Pi^’^ 

=  pu^ . 

{pulni  +  pn^n2)ni, 
{pu^ni  +  /0n^n2)n2, 


Then,  = 

For  the  state  derivatives,  the  symmetry  plane  implies  that  the  normal  derivatives  of 
scalar  quantities — i.e.  p,  pE,  and  pv — are  zero.  The  normal  derivative  of  the  tangential 
velocity  is  zero  as  well.  The  condition  gives  no  information  about  the  normal  derivative 
of  the  normal  velocity  or  the  tangential  derivative  of  any  quantity.  The  boundary  state 
gradient  is  set  using  the  given  information  about  the  normal  derivatives  combined  with 
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interior  derivatives.  Thus,  for  the  density, 


dxi 

dp 

-  dx2  - 


n  -1  r 


ni  n2 
-n2  ni 

ni  n2 
-n2  ni 

ni  n2 
-n2  ni 


-1 


dp 
dn 
dp 

L  ^  J 
0 


-1  r 


^i(^)  +^2ni 
-nms  + 


\9a:2  ) 
dp 
dx2 


nl 


-nin2 


-nin2 


nr 


■ 

dp 

dxi 

dp 

- 

.  dx2  - 

+  ni 


\  dx2  J 


The  energy  and  turbulence  model  working  variable  gradients  are  set  analogously.  The 
velocity  gradients  are  specified  using  =  0  and  interior  gradient  information.  Thus, 
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